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Batch  cation  exchange  and  column  experiments  were 

conducted  to  evaluate  several  exchange  selectivities  which 

have  been  suggested  for  describing  cation  exchange  reactions 

in  solute  transport  models.  The  batch  cation  experiments 

measured  cation  equilibria  with  an  analytical  grade 

macroporous  sulphonic  polystyrene-divinylbenzene  cation 

resin.  The  resin  was  equilibrated  with  binary  and  ternary 

solutions  of  Ca(C104)2-4H20,  Mg (CIO4) 2 • 6H2O,   and  NaC104-H20  at 

constant  CIO4  concentrations.  Adjusted  Vanselow  selectivity 

coefficients  were  calculated  for  exchange  equilibria  with  the 

cation  resin,   as  well  as  for  exchange  equilibria  reported  in 

the  literature  with  montmorillonite  and  with  a  Yolo  loam 

soil.  Least  squares  regression  of  the  natural  logarithm  of 

the  adjusted  Vanselow  selectivity  coefficients  on  the 
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exchanger  phase  cation  equivalent  fraction  were  calculated 
for  both  the  binary  and  ternary  systems. 

The  batch  exchange  and  column  experimental  data  were 
used  to  evaluate  mole  fraction,  equivalent  fraction,  and 
statistical  thermodynamic  exchanger  phase  reference 
functions.  The  predictive  utilities  of  regressions  of  the 
adjusted  Vanselow  selectivity  coefficient  on  exchanger  phase 
composition  were  evaluated.  Four  tests  were  used  to  evaluate 
these  reference  functions  and  regression  equations:  1)  their 
closeness  to  calculated  exchanger  phase  activities;  2) 
accuracy  of  estimates  of  exchanger  phase  composition  based  on 
their  use;   3)   accuracy  of  the  predicted  cation  concentration 
distributions  with  cumulative  column  effluent  volume;  and  4) 
accuracy  of  the  predicted  first  moment  of  the  cation 
distribution  with  cumulative  effluent  volume.  The  mole 
fraction  was  closest  to  the  calculated  exchanger  phase 
activity  for  both  synthetic  resin  and  montmorillonite 
exchangers.  With  a  montmorillonite  system,  the  mole  fraction 
was  the  most  reliable  estimator  of  exchanger  phase 
composition,  while  the  equivalent  fraction  was  the  most 
reliable  for  the  resin  exchanger.  With  the  resin  exchanger, 
the  statistical  thermodynamic  reference  function  gave  the 
most  reliable  estimate  of  column  effluent  cation 
concentration.  In  a  column  packed  with  the  Yolo  loam  soil, 
the    mole  fraction  gave  the  most  accurate  prediction  of 
column  response.  Use  of  variable  Vanselow  selectivity 
coefficients  gave  more  accurate  predictions  in  both  the  batch 
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and  column  experiments, 
the  use  of  ternary  data 
exchange  equilibria  but 


Compared  to  the  use  of  binary  data 

improved  predictions  of  batch 

not  predictions  of  column  response 


INTRODUCTION 


Many  of  the  solutes  which  move  through  natural  and 
artificial  porous  media  are  cations.  Consequently,  cation 
movement  in  soils  and  porous  media  is  worthy  of  study.  The 
movement  of  cations  in  soils  can  affect  soil  fertility,  soil 
salinity,   soil  sodicity  and  the  leaching  of  some  pesticides. 
The  movement  of  cations  in  porous  media  generally  determines 
the  cation  movement  from  brackish  water  injected  into  deep 
aquifers,  the  response  of  ion  chromatography  columns,   and  the 
performance  of  ion  exchangers  used  to  deionize  saline  and 
municipal  waters. 

Cation  movement  in  soils  and  other  porous  media  is 
governed  by  both  physical  and  chemical  processes.  The  major 
physical  process  is  the  mass  flow  of  the  solution  through 
porous  media  which  is  in  turn  governed  by  the  nature  of  the 
porous  media,  the  amount  of  other  fluids   (e.g.  air,   and,  less 
commonly,  water-immiscible  hydrocarbons)  present  and  the 
potential  energy  gradients  within  the  water  which  determine 
the  rate  and  direction  of  the  fluid  flow. 

The  movement  of  cations  through  soil  is  also  affected  by  the 
chemical  interactions  between  solute  in  solution  and  the 
solids  in  the  soil  matrix.  Since  most  soil  and  aquifer 
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minerals  have  net  negative  electrostatic  charge  due  to 
isomorphic  substitution  in  the  mineral  crystal,   ion  exchange 
is  an  important  cation  sorption  mechanism.   In  an  exchange 
reaction,   ions  in  solution  sorb  on  the  exchanger  surface  to 
satisfy  the  electrostatic  charge  of  the  exchanger.  Ions 
adsorbed  by  ion  exchange  can  be  displaced  by  a  counterions 
able  to  satisfy  the  same  amount  of  electrostatic  charge  in 
the  exchanger.  Accordingly,   ion  exchange  reactions  are 
distinguished  by  their  stoichiometry  and  reversibility 
(Helfferich,   1962)  .  Cations  may  be  removed  from  or  added  to 
the  soil  solution  by  mechanisms  other  than  ion  exchange.  For 
example,  the  cation  may  complex  with  organic  ligands  or  sites 
on  the  surfaces  of  soil  minerals.  Another  mechanism  which  may 
control  the  soil  solution  cation  concentrations  is  the 
combined  solubilities  of  the  minerals  present  in  the  soil 
matrix.  For  example,   if  present,   calcium  carbonate  (CaCOs) 
and  other  calcic  minerals  will  control,   to  a  large  extent, 
the  soil  solution  concentration  of  calcium  (Ca2+) .  Another 
mechanism  affecting  soil  monovalent  cations  in  solution  such 
as  potassium   (k1+)   and  ammonium  (NH^)   is  their  formation  of 
inner-sphere  complexes  between  the  layers  of  micaceous 
clays,   again,   if  found,   in  the  soil  matrix.  Once  such 
complexes  are  formed,     these  cations  can  be  released  slowly 
into  the  soil  solution. 


LITERATURE  REVIEW 


Mathematical  Formni  at- i  nn  of  Cat- inn  Tr;^n.c;pr.rl- 


General  partial  differential  equations  which  describe 
the  one-dimensional  convective-dispersive  (sometimes 
referred  to  as  advective-dispersive)  transport  of  each 
cationic  solute  species  i  through  porous  media  under 
transient  water  flow  conditions  may  be  written  in  the  form 
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—   iO  D 


[1] 


where 


z 


-  length  in  the  z  direction   (dm) , 


e 


-  volumetric  water  content   (dm^  dm-3)  , 


D 


=  dispersion  coefficient  (dm2  s"!) , 
=  interstitial  velocity   (dm  s'^) , 


u 


s 


-  ion  i  solution  equivalent  concentration 
(~mol  dm"3)  , 

2  7 


3 
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Pb       =  porous  medium  bulk  density   (kg  dm~3)  , 


e 

m 


,      =  ion  i  solid  solution  equivalent  molality 
(■7-  mol  kg-l)  , 

t         =  time   (s) , 

Lg        =  first-order  rate  coefficient  of  decay  in 
solution   (s~l)  , 

1 

Le        -  first-order  rate  coefficient  of  decay  in  solid 
(s-i) , 

^0 

Lg        -  zero-order  rate  coefficient  of  production  in 


solution   (mol  dm~3  s~^)  ,  and 


0 


Lg        -  zero-order  rate  coefficient  of  production  in 
solid  (s~l) 

(van  Genuchten  and  Alves,    1982)  . 

Since  chemical  processes  such  as  those  described  in  the 
opening  paragraphs  of  the  introduction  may  confound  one 
another  in  both  batch  and  column  experiments,   this  study  will 
restrict  itself  to  the  effects  of  the  ion  exchange  reaction 
upon  the  flux  of  cations  with  a  solution  flowing  through 
porous  media.  To  make  the  examination  of  the  effects  of  ion 
exchange  on  cation  transport  through  porous  media,   the  above 
equation  can  be  simplified  by  restricting  the  system  both 
chemically  and  physically.   Simplification  of  the  equation 
must  be  accompanied  by  concomitant  restrictions  in  the  nature 
of  the  porous  media,   the  chemical  make-up  of  the  solution 
which  flows  through  it,   and  the  physical  nature  of  the 
solution  flow  through  the  porous  medium.  First,  the 
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volumetric  water  content  and  the  dispersion  coefficient  are 
held  constant.  Second,   the  chemical  solute-solid  interactions 
are  limited  to  ion  exchange.  Third,   it  is  assumed  that  the 
ion  exchange  reactions  are  instantaneous  locally.  Ion 
exchange  is  a  kinetic  process,  but  the  rate-determining  step 
of  the  process  is  not  the  reaction  at  the  exchanger  surface, 
but  the  diffusion  of  the  counter-ions  to  and  from  the 
exchanger  surface   (Helfferich,   1962).  The  kinetics  of  the  ion 
exchange  process  will  be  determined  by  the  diffusive  distance 
and  the  resistance  to  mass  flow  the  ion  experiences  as  it 
moves  to  and  from  the  exchanger  surface.  Accordingly,   it  is 
the  structure,   rather  than  the  surface  chemistry,  of  the 
porous  media  which  determines  effective  rate  of  the  exchange 
reaction.     For  porous  media  with  sufficiently  high  molecular 
diffusion  coefficients  and  sufficiently  low  dispersion 
coefficients,  the  ion  exchange  reaction  may  be  assumed  to  be 
occurring  instantaneously  in  solute  transport  models  (James 
and  Rubin,   197  9) . 

If  the  three  restrictions  described  in  the  preceding 
paragraph  are  applied  the  following  general  form  of  partial 
differential  equation  applies  for  each  cationic  solute  i 


S^i'  K''         52  c/,  6c/. 

sir      ~~  ^  ^  67^  "  "  17-  [2] 

The  simultaneous  solution  of  three  of  these  partial 
differential  equations,   one  each  for  three  heterovalent 
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cations  flowing  through  exchanging  porous  media,  will  occupy 
the  balance  of  this  study. 

Thermodynamic  Description  nf  Ton  F.xr.h^ng^ 

By  neglecting  kinetic  effects,  the  ion  exchange 
phenomena  may  be  described  by  equilibrium  thermodynamics. 
Before  beginning  with  this  description,   some  terms  should  be 
defined.  The  exchanaer-sol  nt -i  on  system  is  composed  of  the 
water,  the  solid  ion  exchanger,  which  is  assumed  here  to 
convey  only  a  constant  negative  electronic  charge,   and  the 
salts,  which  are  assumed  here  to  be  completely  disassociated. 
While  the  components  are  commingled,  the  exchanger-solution 
system  can  be  viewed  as  being  composed  of  two  distinct 
phases,  designated  the  exchanger  phase  and  the  solution 
phase.  The  gxghanger  phase  consists  of  the  solid  ion 
exchanger  and  all  cationic  solutes  sorbed  to  the  solid  ion 
exchanger  as  well  as  water  molecules  of  hydration  for  the 
cations  and  the  solid  exchanger.  The  solution  ph^<.^  consists 
of  everything,   aside  from  the  gases  above  the  mixture,  which 
is  not  ascribed  to  the  exchanger  phase   (i.e.  the  water,  the 
anionic  solutes,   and  the  remainder  of  the  cationic  solutes) . 

A  binary  cation  exchange  reaction   (one  which  involves 
two  cationic  solutes)  may  be  expressed  by 


ia:.  .  z,  jx^         z,.  ix:^  .     ja!  [3] 
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where  the  two  cations  i  and  j,   having  valences  Zi  and  Zj, 
respectively,   are  both  assumed  to  be  found  in  the  solution 
completely  disassociated  from  the  anion  Si    and  in  the 
exchanger  phase  sorbed  to  the  exchanger  X.   The  cationic 
solutes  in  the  exchanger  phase  are  indicated  with  a 
superscript  "e"  to  differentiate  them  from  non-exchangeable 
molecules  of  the  same  cation  which  may  be  constituents  of  the 
cation  exchanger. 

The  thgrrnQdynami  r  exchange  constant  for  this  reaction 
relates  the  activities  of  cations  in  the  exchanger  and 
solution  phases: 

Ke^i/j  =    (— )     ^  (~^)  [4] 
a,  a. 

where 

^exi/j  =  thermodynamic  exchange  constant, 

e 

=  activity  of  solute  i  in  the  exchanger  phase, 

and 

3 

=  activity  of  solute  i  in  the  solution  phase. 


Solution  Pha.'.e  and  Exchanger  Pha.se  R^-Ferenr^  Fnno^jr^n^ 

the  activities  of  electrolytes  in  solution  are 

solution  phase  and  exchanger  phase  activities  of 
cationic  solutes  are  not  measurable  directly, 
experimentally  measurable  activity  estimates  are 


While 
measurable, 
individual 
Therefore, 
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used  in  the  application  of  thermodynamic  theory  to 
ion-exchange  reactions  in  soils.  A  reference  function  is  a 
function  of  experimentally  measurable  quantities  which 
approximates  thermodynamic  activities   (Wall,   1974) .  a 
reference  function  must  approach  smoothly  the  activity  as  the 
system  approaches  its  reference  state.  This  criterion  is 
represented  symbolically  by 

lim  aj 

where 

ai       =  activity  of  i, 

r(C)    =  a  reference  function,  and 

C*        =  the  reference  state. 

Both  solution  concentration  and  solution  molality  may  be 
used  as  reference  functions  of  solution  phase  solute 
activity.   Solution  concentration  is  the  appropriate  reference 
function  of  the  solute  in  the  solution  phase  for  models  of 
ion  exchange  with  transport  in  porous  media  when  water 
content  is  expressed  as  a  volumetric  ratio.   The  American 
Society  of  Agronomy   (1984)   has  designated  the  mass  ratio  as 
the  preferred  unit  to  express  soil  water  content.   If  the  mass 
ratio  is  adopted  to  express  soil  water  content,   then  solution 
molality  would  be  a  more  tractable  selection  as  the  reference 
function  for  solution  phase  activity  of  the  solutes. 

Several  exchanger  phase  reference  functions  have  been 
proposed  for  thermodynamic  models  of  ion  exchange.  For 
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purposes  of  this  discussion,   all  exchanger  phase  reference 
functions  can  be  expressed  symbolically  by 

e 

-^i  ~  Di  [5] 

where 

e 

=  exchanger  phase  cation  i  reference  function, 

e 

n_.       =  number  of  cation  i  in  the  exchanger  phase   (mol)  , 

and, 

Di       =  a  denominator,   exact  form  depends  on  particular 
exchanger  phase  reference  function. 

In  1932,  Vanselow  proposed  the  exchanger  phase  mole 
fraction  as  an  estimate  of  the  exchanger  phase  activity.  If 
the  exchanger  phase  activity  is  equal  to  the  exchanger  phase 
mole  fraction   (that  is,   the  exchanger  phase  activity 
coefficient,   fi,   is  unity),   the  exchanger  phase  behaves  as  an 
ideal  solid  solution  obeying  Raoult's  law.  Sposito  and 
Mattigod   (1979)   have  noted  that  trace  metal-Na  exchange  on 
Camp  Berteau  montmorillonite  behave  as  an  ideal  solid 
solution.   It  is  not  clear  if  the  exchanger  phase  of  soils  and 
soil  minerals  behave  consistently  as  ideal  solid  solutions. 
Krishnamoorthy  and  Overstreet   (1948),   however,  have  argued 
that  the  application  of  Raoult's  law  to  ion-exchange 
equilibria  ignores  the  spatial  arrangement  of  charge  sites  on 
the  exchanger  surface.   If  the  mole  fraction  is  accepted  as 
the  exchanger  phase  reference  function,   the  denominator,  D^, 
is 
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n 

e 


Di    -    7  [6] 
k=l 

where 

n       =  number  of  cations  in  the  exchange  system. 


Using  techniques  of  statistical  thermodynamics, 
Krishnamoorthy  et  al.    (194  8)   adopted  the  gas  adsorption 
theory  of  Guggenheim  to  cation  exchange  reactions.  They 
contended  that  the  mole  fraction  reference  function  failed  to 
take  into  account  the  spatial  distribution  of  surface 
functional  groups  on  the  exchanger  surface.   This  shortcoming 
could  be  remedied  by  viewing  the  exchanger  surface  as  a 
two-dimensional  array,   each  surface  functional  group  having 
four  nearest  neighbors.  Krishnamoorthy  and  Overstreet  (1949) 
wrote  that  this  reference  function  would  be  an  appropriate 
model  for  cations  adsorbed  to  "clays  and  synthetic  resins." 
The  result  of  this  approach  suggested  a  more  complicated 
estimate  of  exchanger  phase  activity  with  a  denominator 
appropriate  for  their  estimate  as 


k;=l 


e  +  1 

k    where  yi  =   .  [7] 


Gaines  and  Thomas   (1953)   suggested  the  equivalent 
fraction  as  an  estimate  of  exchanger  phase  activity.  The 
denominator  corresponding  to  this  reference  function  is 
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[8] 


While  the  present  study  ignores  reactions  other  than  ion 
exchange,   there  are  several  computer  programs  in  the 
published  literature  which  attempt  to  calculate  the  complex 
chemical  equilibria  in  soil  solutions.   The  most  sophisticated 
of  these,   GEOCHEM,   is  an  adaptation  of  REDEQL2  for  soil 
solutions    (Sposito  and  Mattigod,    1980;  Morel  and  Morgan, 
1972) .     While  no  explicit  exchanger  phase  reference  function 
is  presented,   the  reference  function  used  implicitly  by 
GEOCHEM  is  the  "concentration"  of  cation-exchanger  complex  in 
solution  in  moles  cation-exchanger  complex  per  cubic 
decimeter  of  solution. 


coefficients.  For  the  purposes  of  this  study,   the  defining 
relation  is  the  equation  which  defines  the  thermodynamic 
exchange  constant  in  terms  of  the  exchanger  phase  and 
solution  phase  reference  functions 


Exchange  .Selertivity  Gnpfficients 


The  exchanger  phase  and  solution  phase  reference 
functions  reviewed  above  may  be  combined  to  yield  selectivity 


■exi/j  - 


[9] 
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where 

^exi/j  =  thermodynamic  exchange  constant, 
e 

=  cationic  solute  i  exchanger  phase  reference 
function, 

s 

=  cationic  solute  i  solution  phase  reference 
function, 

fi        =  exchanger  phase  activity  coefficient,  and 
Yi        =  solution  phase  activity  coefficient. 
The  units  of  the  reference  function  and  the  activity 
coefficient  vary  according  to  the  reference  function 
selected.   Since  activity  is  unitless,  since 


=  fir^  ,  and 


e  s 


the  units  of  the  activity  coefficient  are  the  inverse  of 
those  for  the  reference  function  for  which  it  is  defined. 

The  Inn  exchange  selectivity  coeff  i  mVm-,  ki/j, 
describing  the  exchange  reaction  equilibrium  between  ion  i 
and  ion  j  is  defined  as 


s    /        V   e   /  •  [10 


The  most  recent  recommendation  by  the  International  Union  of 
Pure  and  Applied  Chemistry   (lUPAC)    is  that  the  term 
selectivity  coefficipnt  be  used  when  the  solutes  in  both 
phases  were  expressed  in  terms  of  concentrations  (Irving, 
1972)  .  Since  the  choice  of  the  concentration  scale  was  left 
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up  to  the  investigator  by  lUPAC,  the  most  consistent  usage 
for  this  discussion  is  to  use  the  selected  reference  function 
as  the  appropriate  concentration  scale.  The  adjusted 
selectivity  coeffic-ient ,   k^^.  ,   refers  to  the  selectivity 

coefficient  in  which  the  solution  phase  reference  functions 
are  replaced  by  their  activities 


Most  of  the  commonly  referenced  selectivity 
coefficients,   the  Vanselow  selectivity  coefficient,  the 
Gaines-Thomas  selectivity  coefficient  and  others  are  defined 
with  the  cations  in  the  solution  phase  represented  by  their 
activities  rather  than  a  reference  function.  To  be  consistent 
with  the  recommendations  of  lUPAC,   these  coefficients  will  be 
referred  to  as  adjusted  selectivity  coefficients.  That  is, 
the  common  usage  of  the  selectivity  coefficient  suggested  by 
Vanselow  is  represented  by 


e  s 

^(V)i/j   -   V  s  /      \~i  [12: 

a, 

where 


e 

x_£       -  exchanger  phase  mole  fraction, 
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In  this  study,  this  will  be  referred  to  as  the  "adjusted 
Vanselow  selectivity  coefficient"  )to  be  consistent 

with  the  lUPAC  recommendation,   even  though  this  selectivity 
is  not  different  from  what  is  commonly  understood  to  be  the 
Vanselow  selectivity  coefficient.  Accordingly,   the  term 
"Vanselow  selectivity  coefficient"  )    is  taken  to  be 


e  s 


k 


'i/j       ^    3    '         ^    e  ' 


(V)i/j       ^    s    '         ^    e    ^  [13] 


The  thermodynamic  exchange  constant  is  related  to  the 
selectivity  coefficient  by 


fi  ^3 

and  to  the  adjusted  selectivity  coefficient  by 


^e^^/^-  -       ^zi    ^ilj   ■  [15] 

Other  adjusted  selectivity  coefficients  to  be  considered 

here  include  that  suggested  by  Gaines  and  Thomas  (1953), 
a 

'  the  adjusted  statistical  thermodynamic  proposed 

by  Krishnamoorthy  et  al.(1948),   ;c,3^,_^/^  .  These  adjusted 
selectivity  coefficients  are,  respectively. 


e  s 

^(GT)i/j    -      V    s    /        V~7~"/  [16' 

and 
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C                                    3  n 

.Yi             z-  zi         "V^  z'-z- 

^(ST)i/j    =    \—  )    ^(  —)  ^    (   \yk  ml)  [11] 

yj  k  ^1 


where 

e  ,      ,         ,  , 

X-li  =  cationic  solute  i  exchanger  phase  equivalent 

fraction,  and 


All  of  the  preceding  selectivity  coefficients  come  from 
the  chemistry  of  ion  exchange.  Other  selectivity 
coefficients,   to  be  considered  later,   have  been  suggested  by 
their  inclusion  in  various  solute  transport  models. 

Calculation  of  the  Thermodynamic  Exchange  Constant  and 
Exchanger  Phase  Activity  Coefficients 


The  adjusted  Vanselow  selectivity  coefficients 
determined  from  binary  exchange  experiments  may  be  used  to 
calculate  the  exchanger  phase  activity  coefficients  and  the 
thermodynamic  exchange  constant    (Argersinger  et  al . ,  1950; 
Sposito,    1981a) .     The  thermodynamic  exchange  constant  is 
related  to  the  adjusted  Vanselow  selectivity  coefficient  by 

f. 

=       \  Zi  [18] 
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which  by  taking  the  natural  logarithm  of  both  sides  yields 

In  i^exi/j  =  In  ^(vU/j  +         1"        -         ^^^j-  [19] 

At  equilibrium,   changes  in  chemical  potential  of  two 
cations  and  water  in  the  exchanger  phase  are  related  by  the 
appropriate  form  of  the  Gibbs-Duhem  equation 


where 
e 

-  number  of  moles  of  cation  i  in  the  exchanger  phase, 

e 

=  number  of  moles  of  j  in  the  exchanger  phase, 

=  number  of  moles  of  solvent  in  the  exchanger  phase, 

e 

=  exchanger  phase  cation  i  chemical  potential, 
e  _ 

\ij       -  exchanger  phase  cation  j  chemical  potential,  and 

e 

^„       =  exchanger  phase  solvent  chemical  potential. 


While  it  is  not  essential  to  the  validity  of  the 
resulting  equations,    it  will  be  assumed  here  that  d^J  =  0 .  By 
using  the  definition  of  the  chemical  potential  of  the  two 
cations  in  the  exchanger  phase 
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e        eo  e 
^i^  =         +  RT  In  x^fi,  and  [21] 

e        eo  e 
[Lj  =  [L  J    +  RT  In  xjfj, 

where 
eo 

[1  ^       =  standard  state  exchanger  phase  cation  i  chemical 
potential,  and 

eo 

|J.  J       =  standard  state  exchanger  phase  cation  j  chemical 
potential . 

The  standard  state  for  this  discussion  will  be  homoionic 
exchanger  in  contact  with  a  solution  having  the  same  anion 
concentration  as  the  equilibrating  solutions  and  only  cation 
i  or  j.  By  using  the  fact  that, 

dxi  =  -dxj 
the  following  equation  may  be  derived: 

Xi'dlnfi  +  Xj 'din  fj    =0  [22] 
By  using  the  derivative  of  the  relation 

In  iCexi/j  =  In  k^^^.^^  +  zj  In  fi  -  zi  In  fj 
which  is 


'^^^^  )   =        dilnfj)    -  Zj  d(ln  fi  ),  [23] 


the     differential  equations 

Zj  dlnfi  =   (1-x/,  )x^.  dln;c,J,       ,   and  [24] 
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zidlnfj  =  X..  ciln;c,J,,/.  [25] 


where 


e 

Xj^,  -  exchanger  phase  cation  i  equivalent  fraction,  and 
e 

Xj,  -  exchanger  phase  cation  j  equivalent  fraction 
may  be  derived.   These  equations  may  be  integrated  to  yield 

1 


and. 


in  fj  =  jj  (x/,  ln;c,J,,/.  -  Jlnk^;^.^.  dx^.  )  [27] 


Recalling  that 


In  Kexi/j  =  In  k^^^.^.  +  z^-  In  fi  -        In  fj 

these  equations  may  be  combined  to  provide  an  estimate  of  the 
thermodynamic  exchange  constant 


1 

e 


[28] 

6 

The  standard  free  energy   (AG°)   of  the  exchange  reaction  may 
be  calculated  from  the  exchange  constant  by  the  use  of  the 
relation : 

AG°  =   -RT   In   Kexi/j  , 
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where 


R 


Molar  gas  constant    (J         mol"^)  , 


T 


Temperature   (K)  , 


(Lewis  and  Randall,   1961) .  This  approach  has  been  validated 
for  the  calculation  of  the  thermodynamic  exchange  constant 
for  soils    (See,    for  example,   Jensen  and  Babcock,    1973)  . 

The  Cation  Transport  Model  of  Rubin  and  James 

The  first  numerical  model  of  multispecies  cation 
transport  through  porous  media  to  gain  wide  acceptance  was 
that  of  Rubin  and  James   (1971) .   This  seminal  article  is  both 
a  good  example  of  integrating  complex  chemical  submodels  in 
solute  transport  models  and  an  example  of  choosing  the  form 
over  the  substance  in  the  selection  of  a  chemical  submodel. 
Recalling  the  convect ive-dispersive  (or 

advective-dispersive)   equation   [2]   from  the  introduction 


+ 


=  D 


-  u 


5t 


0  5t 


6t  2 


which  may  be  written  as 


5t 


+ 


e  5t 


[29] 


where  the  operator  L  is  equal  to 
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52  5 

L  =  D  - — r  -  u  — 
5t  2  5^ 


Rubin  and  James  defined  their  exchange  selectivity 
coefficient  as 

e  3 
,  _  {_±  \  ^3  I    0  \ 

^(Rj)i/j      V    3    /        \    e    /         •  [30] 

Although  Rubin  and  James  state  that  this  exchange  selectivity 
coefficient  need  not  be  considered  a  constant,   efforts  to 
include  variable  select ivit ies  in  models  following  their 
development  have  not  proven  altogether  successful  for 
selected  data.  Most  authors    (e.g.  Valocchi  et  al . ,   1981b)  who 
have  followed  development  by  Rubin  and  James  have  assumed 
constant  exchange  selectivity. 

In  their  development,   Rubin  and  James  made  two 
assumptions.  First,   they,   as  do  most  authors,   assumed  that 
the  cation  exchange  capacity  (C  )   of  the  porous  medium  and 
the  selectivity  coefficients    (k^^J,ifj  )   for  each  pair  of 
cation  species  are  constants. 


i  -1 


They  define  the  function,  Uij,  a  multispecies  exchange 
isotherm  as 


^7  ^  i  z  .  z  v 

ij      K  (Rj)  ij  c.,        rn  ^ ,      _  m  e         ^  s 


[31] 
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which  after  partial  differentiation  with  respect  to  time,  t, 

5ui  -j 

and  division  by  — ^  yields 

6m  j , 


m  ^ 


dt  du  at 

 £j 

dm^  ,  a  771  J  , 


9cs  a::^  dn  ^ , 

r  r  J 

+      +    =  0 

du  dt  dt 


noting  that 


1  — 

j  =  1     3  t 


d  m  j 


—  0 


and  that 

dm.,  -    .  . 


d  m  J  ,  d  m  f 


dt  dt 


[32] 


and  summing  the  partial  different  equation  above  for  all 
j=l,n  gives 
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d  m 


3t  '         at  at 


where 


and 


du 

 ij 

j  *i     ou  .  , 

 IJ 

d  m  ^ , 


du 

-  i  ' 


 £j 

d  m  ^ , 


r 


V  .  . 


du 


du 


8  /n  ^ , 


[33] 


which,   when  substituted  into  the  equation   [29],   gives  the 
relation 

f  o  \  d  c  \  '  -r-.  d  c\' 

(9  -  P   P  y  ^ij   —    =  L  c;, 

at  ^  dt 

[34] 
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This  equation  is  then  solved  numerically,   under  the 
appropriate  boundary  and  initial  conditions,  to  give  changes 
in  cation  solution  concentration  with  time. 


Crank-Nicholson  Method 


For  the  one-dimensional  flow  regime,  equation   [34]   is  a 
parabolic  partial  differential  equation  which  may  be  solved 
numerically  with  a  finite  difference  method,  that  is  the 
Crank-Nicholson.  Finite  difference  methods  are  relaxation 
methods,  which  means  that  the  errors,  or  residuals,  resulting 
for  an  initial  approximation  are  considered  as  constraints 
that  are  to  be  relaxed  (Hornbeck,   1975) .  New  approximations 
are  chosen  to  reduce  the  worst  of  the  residuals  until  finally 
all  are  within  a  tolerance  limit.  The  finite  difference 
methods  are  based  on  the  Taylor  series  expansion  where  the 
value  of  the  function  u  at  the  point  j  +  Ax  is  approximated 
by 


u  ( j  +  A  X  )  =    u  ( j  )  +  A  X  u   '  ( j  )  +  -^-li  u   "  {j  )  + 

2! 

(James  and  James,   1976) . 

The  first  derivative  of  the  function  u  with  respect 
to  X  is  therefore 

u     ij  ;  -   ■    u  "ij  )  + 


Ax 


2! 
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which  may  be  expressed  as 

u  ( j  +  A  X  )  -  u  ( j  ) 


u  '  ( j  )  = 


Ax  Ax 

where 

R2  =  is  the  remainder  after  2  terms  in  the  Taylor 
series . 

Similarly,   the  second  derivative  may  be  estimated  by 


„  ,  .  ,  u  (J  +  A  X  )  -  2  u  I  j  )  +  u  (j  -  A  X  )  R3 
u  "  {J  )  =  1  I   -  


Ax2  Ax2 

There  are  several  possible  finite  difference  solutions 
of  parabolic  partial  differential  equations;  the 
Crank-Nicholson  is  frequently  used  because  the  estimate  is 
second-order  accurate. 

For  example,   the  numerical  approximation  of  the 
parabolic  partial  differential  equation 


dx2  By 

[35] 

would  be  expressed  in  terms  of  the  values  of  the  five  nearest 
neighbors  of  the  node  of  interest.  The  parabolic  equatic 
[35]   would  be  approximated  by  the  equatic 


-on 

Lon 
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_a_  r  f  {j  +  ^x,k+Ay)_2f(j+Ax,k)^f  {j  +  Ax  ,k  -  Ay  ) 
2  (  Ax)^ 


+    -f  (J  '  k  +  Ay  )  -  2  f  {j  ,  k  )  +  f  {j  ,  k  -  Ay  ) 


(  Ax 


f   ij  +  A  X ,  k  >    -  f   (j  ,  k  ) 


(  Ay  )  ' 


This  equation  for  the  x-y  region  of  interest  can  be 
expressed  in  matrix  form.   The  resulting  matrix,   in  turn,  can 
be  solved  in  a  straightforward  manner  using  matrix 
multiplication . 


1    r    1  f  *  2 


a 


1     Y  1 


j  +1,  2 


j  +1,  3 


0 


j  +1, n  +1  ^  n  -1 


where 


] 
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a  (  Ay  ) 


L        a  (  Ay  )    J   ^  ' 


<I>     -    I  <i  -    I  _f         -    f  -  f 

L  a(Ay)J^'*  j  .  k  -1      ^  j  .  k  +1 


Exchanger  Phase  Peference  Function.^  Tmpiiri^  jn  p^m^ 

Transport  Modpl.^ 


Earlier  in  this  literature  review,  exchange 
selectivities  and  exchanger  phase  reference  functions  which 
have  been  developed  for  the  chemical  study  of  exchange 
reactions  were  discussed.   This  section  reviews  the  chemical 
submodels  of  solute  transport  models  which  simulate  the 
transport  of  cationic  solutes  through  reactive  porous  media. 
A  chemical  submodel  of  a  solute  transport  model  is  defined 
here  as  that  part  of  the  solute  transport  model  which 
describes  the  chemical  reactions  involving  the  solute.  Most, 
though  not  all,   of  these  chemical  submodels  of  transport 
equations  incorporate  mass  action  equations  and  exchange 
selectivities.  The  exchanger  phase  reference  functions  which 
are  implicit  in  the  exchange  selectivities     used  in  these 
chemical  submodels  are  typically  different  from  those  current 
in  chemical  reports  of  exchange  reactions.  Though  all  types 
of  cationic  solute  chemical  submodels  will  be  reviewed  here, 
most  of  the  attention  will  be  placed  on  chemical  submodels 
which  define  exchange  selectivities. 
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The  chemical  submodels  found  in  published  models  of 
cation  transport  through  porous  media  can  be  separated  into 
three  general  categories:  linearized,   surface  complexation, 
and  exchange  selectivity  coefficient.  Many  sorption  phenomena 
in  soils  have  been  found  to  approximate  linearity  and  can  be 
described  successfully  by  a  partition  coefficient,  K^. 
Several  authors    (Nkedi-Kizza  et  al . ,    1984;  Valocchi,  1984; 
Miller  and  Benson,   1983),   while  noting  that  ion  exchange  is 
freguently  a  non-linear  phenomenon,   have  suggested  that 
under  certain  conditions  ion  exchange  reactions  in  solute 
transport  models  could  be  described  by 


s  e 

[36] 


The  partition  coefficient  could  be  calculated  from  batch 
experiments   (Nkedi-Kizza  et  al . ,    1984)   or  estimated  from  the 
exchange  selectivity  coefficients  based  implicitly  on  either 
mole  fraction   (Miller  and  Benson,   1983)   or  equivalent 
fraction   (Valocchi,   1984)   reference  functions. 

Recently,   surface  complexation  models  have  been  applied 
to  models  of  cation  transport.  Cederberg  et  al .  (1985) 
incorporated  the  surface  complexation  models  available  in 
MICROQL  to  solute  transport  equation  named  TRANQL .  Jennings 
et  al.    (1982)  applied  the  surface  complexation  model  of 
Schindler  et  al .    (1976)   in  their  cation  transport  model. 
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The  thermodynamic  exchange  constant  relates  the 
activities  of  the  two  cations  in  the  exchanger  and  solution 


phases 


a,  a. 


where 


^exi/j  =  thermodynamic  exchange  constant 


e 


-  activity  of  solute  i  in  the  exchanger  phase 
s 

=  activity  of  solute  i  in  the  solution  phase 

=  valence  of  cationic  solute  i 
By  using  the  mathematical  form  of  the  thermodynamic 
exchange  constant  defining  relation  in  establishing  a 
chemical  submodel,   the  developers  of  cation  transport  models 
are,   consciously  or  unconsciously,   using  the  resulting 
implicit  exchanger  phase  reference  function  as  an 
approximation  of  exchanger  phase  activity.  Mole  fraction, 
equivalent  fraction,   equivalent  molality,   Gapon,   and  mass 
ratio  have  all  been  used  as  exchanger  phase  reference 
functions . 

Many  authors  have  followed  the  seminal  work  of  Rubin  and 
James   (1971)   who  defined  their  exchange  selectivity 
coefficient    (using  the  symbols  of  this  study)  as 
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This  coefficient  uses  the  equivalent  exchanger  phase 

molality   (j-  mol  kg-i)   as  the  implicit  exchanger  phase 

reference  function.  This  implicit  reference  function  has  been 
adopted  also  by  Robbins  et  al.    (1980),   Cederberg  et  al . 
(1985),   Charbeneau   (1981),   Barnes  and  Aylmore   (1984),  Rubin 
(1983),   and  Dutt  et  al .    (1972).  Among  the  selectivity 
coefficients  used  in  soil  chemistry,  the  Gaines-Thomas 
reference  function  is  closest  to  the  selectivity  coefficient 
used  by  Rubin  and  James.  The  two  selectivity  coefficients  are 
related  by  the  relationship 


■(GT)i/j  =    Zi'^^Zj  ;C(Rj)i/j   CEC  ^ 


[37] 


where  CEC  is  the  cation  exchange  capacity  of  the  porous 
medium. 

Some  authors  have  chosen  the  equivalent  fraction  as  the 
exchanger  phase  reference  function  for  exchange  chemical 
submodels  of  solute  transport  models   (Valocchi  et  al.  1981a 
and  1981b;  Lai  and  Jurinak,   1971;  and  Meyers  and  Salter, 
1984) . 

If  we  denote  the  exchange  selectivity  coefficient 
defined  by  the  choice  of  the  exchanger  phase  equivalent 
fraction  and  the  solution  phase  normality  as  the  Valocchi 
selectivity  coefficient  (;C(VAL)i/j) 
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it  is  easy  to  show  that  this  selectivity  coefficient  is 
related  to  the  Rubin-James  selectivity  coefficient  by 

^(VAL)  i/j  =  -^(Rj)  i/j  CEC^^'^i .  [39] 

Accordingly,  the  exchanger  phase  equivalent  molality  and 
exchanger  phase  equivalent  fraction  will  be  treated  in  this 
study  as  being  functionally  equivalent.  For  reasons  of 
brevity,   only  the  exchanger  phase  equivalent  fraction  and  the 
Gaines-Thomas  selectivity  coefficient  will  be  discussed, 
while  the  conclusions  drawn  for  these  representations  apply 
directly  to  representations  based  on  the  exchanger  phase 
equivalent  molality  reference  function. 

While  the  mole  fraction  reference  function  and  the 
Vanselow  selectivity  coefficient  which  it  defines  have  been 
used  widely  in  investigations  of  ion  exchange  equilibria  in 
batch  systems,   only  one  cation  transport  model   (that  used  by 
Miller  and  Benson,   1983)   explicitly  uses  that  reference 
function . 

Apparently  most  developers  of  cation  transport  models 
have  found  that  the  exchanger  phase  equivalent  fraction 
reference  function   (or  its  analogue)   is  more  straightforward 
to  use  in  cation  transport  models  than  the  exchanger  phase 
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mo. 


)le  fraction  reference  function.  This  occurs  because  if  the 
cation  exchange  capacity  of  the  porous  medium  is  assumed  to 
be  constant  and  the  quantity  of  exchangeable  cations 
expressed  in  numbers  of  equivalents,  the  sum  of  exchangeable 
.  cations  per  unit  mass  of  the  porous  medium  will  always  equal 
to  the  CEC,  whereas  the  total  number  of  moles  of  exchangeable 
cations  will  vary  with  exchanger  phase  composition.  This 
difficulty  could  be  overcome  if  the  Gaines-Thomas 
selectivity  coefficient   (or  its  analogue)  were  allowed  to 
vary  as  a  function  of  exchanger  phase  composition.  Even 
without  considering  variable  Gaines-Thomas  selectivity 
coefficients,  these  solute  transport  models  are  often  very 
complex,   so  that  the  apparent  inclination  has  been  to  hold 
the  Gaines-Thomas  selectivity  coefficient   (or  its  analogue) 
constant.  Only  Lai  an  his  coworkers   (Lai  and  Jurinak,  1972; 
and  Lai  et  al . ,   1978)   and  Mansell  et  al.    (1986)  have 
considered  exchange  selectivities  which  vary  with  exchanger 
phase  composition. 

Thus  two  questions  must  be  considered  when  models  are 
used  to  describe  cation  transport  in  porous  media:  1)  what  is 
the  magnitude  of  error  conveyed  to  predictions  of  cation 
movement  in  porous  media  by  adoption  of  a  particular 
exchanger  phase  reference  function  and    2)     to  what  extent 
are  transport  model  simulations  improved  by  using  exchange 
selectivities  which  vary  with  exchanger  phase  composition. 
This  investigation  was  designed  to  provide  answers  for  these 
questions  using  a  simplified  system  of  three  cationic  solutes 
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in  coluitms  packed  with  a  mixture  of  quartz  sand  and  synthetic, 
exchange  resin. 


OBJECTIVES 


It  was  assumed  that  all  exchanger  phase  reference 
functions  could  describe  exchange  behavior  equally  well  if  a 
functional  relation  could  be  developed  for  the  chosen 
reference  function  and  exchanger  phase  composition.   If  the 
selectivity  coefficient  was  set  to  a  fixed  value,  however, 
which  selectivity  coefficient  would  give  the  most  accurate 
description  of  exchange  equilibria?  The  primary  objective  of 
this  study  was  to  determine  which  exchanger  phase  reference 
functions,  when  employed  in  a  chemical  submodel  of  a  solute 
transport  model,  gave  the  most  accurate  description  of  cation 
movement  in  columns  of  porous  media  in  systems  involving 
ternary  cation  exchange. 

The  secondary  objective  was  to  determine  if  data  from 
binary  ion  exchange  experiments  were  sufficient  for  the 
accurate  description  of  heterovalent  ternary  exchange 
isotherms  and  the  accurate  description  of  cation  movement  in 
columns  of  porous  media  in  systems  involving  ternary  cation 
exchange . 

Three  exchanger  phase  reference  functions  were  evaluated 
in  this  study:  the  mole  fraction,  the  equivalent  fraction, 
and  the  statistical  thermodynamic.  The  predictive  utility  of 
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these  reference  functions  were  to  be  evaluated  for  three 
capabilities:  1)   to  estimate  accurately  exchanger  phase 
activity  in  binary  exchange  experiments;  2)  when  used  to 
calculate  binary  and  ternary  cation  exchange  equilibria,  to 
estimate  accurately  exchanger  phase  compositions;  and  3)  when 
incorporated  in  the  chemical  submodel  of  a  solute  transport 
model  simulating  the  one-dimensional  flow  of  three  cations 
through  a  column  packed  with  exchanging  porous  media,  to 
predict  accurately  the  distribution  of  cations  in  the  column 
effluent . 

Typically,  binary  exchange  experiments  must  be  conducted 
to  determine  the  thermodynamic  exchange  constant.  For 
accurate  modeling  of  a  ternary  or  higher  order  exchange 
system,  however,  the  exchange  constant  alone  is  not 
sufficient  and  multicomponent  exchange  experiments  should  be 
conducted  to  augment  the  information  obtained  by  binary 
exchange  experiments.  An  addition  objective  of  this  study  was 
to  assess  the  increase  in  predictive  accuracy  ternary 
exchange  information  brought  when  used  to  calculate  exchange 
equilibria  both  in  batch  ternary  exchange  experiments  and  in 
simulations  solute  concentrations  in  column  experiments  with 
three  cations. 


MATERIALS  AND  METHODS 


Batch  Determinririon  of  Binary  and  r^m^ry  Exrhanap  Tcnthorni- 

MatRr-i  a1 

The  ion  exchanger  used  for  these  experiments  was  AG  MP- 
50  Analytical  grade  Macroporous  Cation  Resin,   200-400  mesh, 
which  was  supplied  by  BioRad  in  hydrogen  form.  AG  MP-50  is  a 
sulphonic  polystyrene-divinylbenzene    strongly  acidic  cation 
exchanger  with  8%  crosslinking  marketed  by  BioRad.  The  listed 
exchange  capacity  was  4.6  mol(+)   kg-i  with  a  nominal  surface 
of  area  of  35  m2  g-1 .  The  lyotropic  series  of  this  exchanger 
for  monovalent  and  divalent  cations  are   (Dow  Chemical 
Company,   1964) : 

Ba2+>pb2+>Sr2+>Ca2+>Ni2+>cd2+>Cu2+>Co2+>Zn2+>Mg2+>Be2+ 

and 

Agl+>Til+>Csl+>Rbl+>Kl+>NH4l+>Nal+>Hl+>Lil+. 
By  comparison,   the  general  lyotropic  series  of  soils  and  clay 
minerals  for  monovalent  and  divalent  cations  are 

Ba2+>Sr2+>Ca2+Mg2+ 

and 

Hl+>Agl+>Til+>Csl+>Rbl+>Kl+=NH4l+>Nal+=Lil+ 
(Bohn  et  al . ,    1985) . 


35 


36 


The  gross  and  fine  microscopic  structure  of  this  type  of 
exchanger  are  presented  in  Figures  1  and  2.  For  this  type  of 
exchanger  at  25°,   the  standard  free  energy  of  the  exchange 
reactions  Ca->Na   (meaning  that  the  reaction  product  is  Na-ion 
exchanger),  Mg-^Ca,   and  Na-^Mg  were  found  to  be  1500,  -2800, 
and  -170  J  mol-l,   respectively   (Boyd,   1970) .  Perchlorate  was 
chosen  as  the  common  anion  because  unlike  chloride, 
perchlorate  does  not  form  monovalent  ion  pairs  or  complexes 
in  solution  with  bivalent  alkaline  earth  cations  such  as 
calcium  and  magnesium  (Sposito  et  al . ,    1983)  .  Three 
perchlorate  salts  were  used  in  the  experiments  calcium 
perchlorate   (Ca  (CIO4) 2 • 4H20,   formula  weight  371.06,  supplied 
by  GFS  Chemicals),  magnesium  perchlorate   (Mg  (CIO4) 2 • 6H2O, 
formula  weight  331.33,   supplied  by  GFS  Chemicals),   and  sodium 
perchlorate   (NaC104-H20,   formula  weight  140.46,   supplied  by 
Fisher  Scientific) . 

Seven  exchange  isotherms,   four  ternary  and  three  binary, 
were  determined  by  batch  techniques.  As  discussed  in  the 
literature  review,  the  value  of  the  thermodynamic  exchange 
constant  may  be  estimated  by  integrating  the  values  of  the 
adjusted  Vanselow  selectivity  coefficient  from  the  two 
homoionic  ion  exchanger  compositions,   the  exchanger  phase 
equivalent  fraction  acting  as  the  variable  of  integration. 
Accordingly,   equilibrating  solution  compositions  were  chosen 
so  that  equilibria  were  established  at  regular  intervals  of 
the  exchanger  phase  equivalent  fraction.   Since  the  exchange 
constants  were  not  known,   the  expected  exchanger  phase 
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Polystyrene  chain 
Divinylbenzene  cross-linking 
Sul phonic  group 


Figure   2.    Schematic  drawing  showing 
the  microscopic  structure  of  a 
sulphonic  divinylbenzene-polystyren. 
cation  resin.    (After  Paterson,  1970) 
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compositions  were  only  approximate.   In  establishing  the 
equilibrating  solution  compositions,   it  was  assumed  that 
sodium  was  the  cation  least  preferred  of  the  three  cations  by 
the  exchange  resin  and  the  ion  exchange  resin  showed  equal 
preference  for  calcium  and  magnesium  ions.  Target 
equilibrating  solution  cation  equivalent  fractions  are 
presented  in  Tables  1  and  2. 

Experimental  Pror.(='d^irf^ 

The  batch  and  column  experiments  were  performed  in  a 
laboratory  in  which  the  ambient  temperature  was  maintained  at 
24.5  ±  0.6°  C.  Each  equilibration  was  replicated  thrice.  Each 
batch  equilibration  was  conducted  in  a  50  cm3  polycarbonate 
Oak  Ridge  centrifuge  tube.  After  the  centrifuge  tubes  were 
oven  dried  and  allowed  to  cool  over  anhydrous  calcium  sulfate 
in  an  evacuated  desiccator,   the  combined  mass  of  each  marked 
tube  and  its  cap  was  determined  to  a  0.1  mg  on  a  Mettler  H80 
balance.  Approximately  0 . 1  g  of  AG  MP-50  cation  resin  was 
placed  in  each  tube.   Since  the  mass  of  the  equilibrating 
solutions,  described  in  following  paragraph,  was 
approximately  20  g,  the  solids  to  solution  ratio  for  these 
experiments  was  roughly  1:200.  The  first  series  of 
equilibrations  was  conducted  to  replace  exchangeable  hydrogen 
on  the  exchange  sites  with  a  cationic  composition  close  to 
that  of  the  final  equilibrations.   The  cation  resin  was 
equilibrated  with  solutions  prepared  with  the  same  rational 
composition  as  indicated  in  Tables  1  and  2,   but  with  a 
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Table  1.  Target  cationic  equivalent  fractions  for  the  ternary 
exchange  equilibrating  solutions. 


 Cation 

Set 


number        Equilibration  Ca 


Mg  Na 


-equivalent  fraction- 


A  0.720  0.180  0.100 

B  0.520  0.130  0.350 

C  0.280  0.070  0.650 

°  0.120  0.030  0.850 

E  0.100  0.025  0.875 

F  0.080  0.020  0.900 

G  0.060  0.015  0  925 

H  0.020  0.005  0  975 


A  0.540  0.360 


0.100 

B  0.390  0.260  0.350 

^  0.210  0.140  0.650 

°  0.090  0.060  0.850 

E  0.075  0.050  0.875 

^  0.060  0.040  0.900 

^  0.045  0.030  0  925 

^  0.015  0.010  0.975 


Three  a  0.3  60 

B  0.260  0.390 

C  0.140  0.210 


0.540  0.100 
0.350 


D  0.060  0.090 

E  0.050  0.075 

F  0.040  0.060 

G  0.030  0.045 

H  0.010  0.015 


Four  A  0.180 

B  0.130  0.'520 

C  0.070  0.280 


D  0.030  0.120 

E  0.025  0.100 

^  0.020  0.080 

G  0.015  0.060 

H  0.005  0.020 


0.  650 
0.850 
0.875 
0.900 
0.925 
0.  975 


0.720  0.100 
0.350 


0.650 
0.850 
0.875 
0.900 
0.925 
0  .  975 
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Table  2.  Target  cationic  equivalent  fractions  for  the  binary 
exchange  equilibrating  solutions. 


Set 

number 


Equilibration 


Ca 


Cation 


Mg 


Na 


Five 


A 
B 
C 
D 
E 
F 
G 
H 
I 
J 
K 


 equivalent  fraction- 

0.000  1.000 

0.100  0.900 

0.200  0.800 

0.300  0.700 

0.400  0.600 

0.500  0.500 

0.600  0.400 

0.700  0.300 

0.800  0.200 

0.900  0.100 

1.000  0.000 


Six  A  0.000 


Seven  a 


000 

B                  0.025  0  975 

C                  0.050  0^950 

D                   0.075  0  925 

E                   0.100  0^900 

F                   0.125  0.875 

G                   0.150  0.850 

"                  0.350  0.650 

I                  0.650  0.350 

J                  0.900  0.100 

K                  1.000  0.000 

0.000  1.000 


S  0.025  0.975 

^  0.050  0.950 

^  0.075  0.925 

^  0.100  0.900 

F  0.125  0.875 

^  0.150  0.850 

H  0.350  0.650 

I  0.650  0.350 

J  0.900  0.100 

^  1.000  0.000 
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constant  chloride  concentration  of  1.0  mol  dm-3.  The  cation 
resin  was  equilibrated  with  40  cm^  of  the  appropriate 
solutions.  Because  the  effective  mean  diameter  of  the  resin 
was  large  as  compared  with  soil  colloids,  the  resin  settled 
to  the  bottom  of  the  centrifuge  tubes  after  a  few  moments, 
and  centrifugation  was  not  necessary.  Once  the  resin  had 
settled  to  the  bottom  of  the  centrifuge  tube,   the  supernatant 
solution  was  decanted  and  a  fresh  20  cm3  of  solution  were 
added  to  the  centrifuge  tubes.  The  equilibration  and 
decantation  were  repeated  until  the  pH  of  the  supernatant 
solution  was  greater  than  6.0,   which  was  achieved  typically 
after  seven  equilibrations.  After  this,   the  pH  of  the 
equilibrating  solution  was  not  monitored  further.  While 
Sposito  and  Fletcher   (1986)   have  argued  that  H+  may  be 
considered  as  an  exchangeable  cation  in  some  soil  systems, 
the  pKa  of  the  sulphonic  surface  functional  group  of  the 
cation  resin  is  2  and  thus  it  may  be  safely  assumed  that  H+ 
cations  were  absent  from  the  exchange  complex  during  the 
column  and  batch  experiments    (Kunin,   1972) .  The  resin  was 
then  equilibrated  with  four  rinses  of  deionized  water.  The 
resin  was  then  equilibrated  eight  times  with  20  cm-3  of  salt 
solutions  having  rational  cationic  compositions  indicated  in 
Tables  1  and  2,   but  with  a  constant  perchlorate  concentratic 
of  0.01  mol  dm-3.   The  supernatant  solution  of  the  final 
equilibration  was  retained  and  the  combined  mass  of  the  cap 
and  the  centrifuge  tube  containing  the  resin  and  entrained 
perchlorate  solution  was  measured  to  an  accuracy  of  0.01  g 


Lon 
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with  a  Mettler  digital  balance.  The  centrifuge  tubes 
containing  the  resin  and  entrained  solution  were  then 
oven-dried,   allowed  to  cool  over  calcium  sulfate  in  an 
evacuated  desiccator  and  weighed  to  the  nearest  0.1  mg.  The 
resin  was  then  extracted  by  five  washings  of  20  cm^  i.o  mol 
dm-3  ammonium  acetate   (NH4C2H3O2)  .  The  supernatant  from  each 
washing  was  poured  into  a  100  cm^  volumetric  flask,   which  was 
brought  to  volume  after  the  final  washing.  The  extracted 
solutions  were  filtered  to  remove  any  resin  which  found  its 
way  into  the  volumetric  flask  and  stored  in  125  cm^  Nalgene 
rectangular  high-density  polyethylene  bottles  until  analysis. 

Equilibrium  and  extracted  solutions  were  diluted  in 
duplicate  to  yield  solutions  to  analyze  for  calcium, 
magnesium,   and  sodium.  All  chemical  analyses  of  the 
equilibrium  and  extractant  solutions  were  performed  on  a 
Perkin-Elmer  4  60  Atomic  Absorption  Spectrophotometer.  Calcium 
in  solution  was  determined  by  atomic  absorption  flame 
spectrometry   (Willard  et  al.,   1981;  Perkin-Elmer  Corporation, 
1976) .  A  Perkin-Elmer  Intensitron™  Lamp   (#303-6017)   was  the 
light  source.  A  nitrous  oxide-acetylene  flame  flame  was  used. 
Absorbance  was  measured  at  422.7  nm.  Magnesium  was  analyzed 
by  atomic  absorption  flame  spectrometry.  An  acetylene-air 
flame  was  used.   A  Perkin-Elmer  Intensitron™  Lamp  (#303-6042) 
was  the  light  source.  Absorbance  was  measured  at  285.2 
Sodium  was  determined  by  flame  emission  spectrometry.  An 
acetylene-air  flame  was  used.  Emission  was  measured  at 
589.6  nm. 


nm, 
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Data  Analygp.c; 

Gravimetric  data  were  used  to  calculate  the  weight  of 
the  oven  dry  cation  resin  and  the  weight  of  the  entrained 
solution.  The  chemical  analyses  yielded  the  concentrations  of 
calcium,  magnesium,   and  sodium  in  the  equilibrating  solutions 
and  the  extractant  solutions.  Concentrations  of  the  ions  in 
the  equilibrating  solutions  were  taken  to  be  that  in 
entrained  solutions.  The  following  formula  was  used  to 
calculate  the  molalities  of  each  ion  in  the  exchanger  phase: 


V{ 


(v.  f . )      (extr. ) -V(entr.)      (equil . ) 


m.  = 


^  m(resin) 
where 

V(v.f.)       =  volume  of  the  volumetric  flask  (dm3) 

V(entr.)      =  volume  of  the  entrained  solution  (dm^) 
s 

c_.(extr.)    =  concentration  of  cation  i  in  the  extractant 
solution   (mol  dm~3) 

s 

(equil.  )=  concentration  of  cation  i  in  the  final 
supernatant  equilibrating  solution 
(mol  dm~3) 
m(resin)      =  mass  of  the  resin 

The  mean  liquid  phase  and  surface  phase  concentrations 
for  the  three  replications  were  used  to  for  all  further 
calculations.  The  liquid  phase  activity  coefficient  is 
estimated  by  the  Davies  equation 


1/2 

Iny,    =  -  A^z/  (   J,  ) 


where 

Yi  =  is  the  solution  phase  cation  activity  coefficient, 
Ay  =  Debye-Hiickel  parameter, 

Ic  =  ionic  strength  of  the  electrolyte  solution, 
(concentration  basis) 

-it  ^  . 

k  =1 

Zjt  =  valence  of  cation  or  anion. 

(Sposito,   1981b) . 
These  data  were  used  to  calculate  Vanselow,  Gaines-Thomas, 
and  statistical  thermodynamic  adjusted  selectivity 
coefficients . 

Quadratric  regressions  were  performed  for  the  relation 
between  In  ,    Inic,,?,,/^  ,   ln;c,s?,i/^    and  the  surface 

phase  equivalent  fraction  of  the  cation  j.  The  regression 
between  Ink  and  the  surface  phase  mole  fraction  was 

used  to  evaluate  the  integral 

1 

In  i^exi/j    =  Jm^r.J,,/.  dx^l  . 

0 

Binary  and  ternary  regression  parameters  were  adjusted 
so  that  the  triangle  rule  appropriate  for  thermodynamic 
exchange  constants  was  obeyed   (Helfferich,  1962). 
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Cglumn  Experiments  Used  to  A.g.^ess  ChPmiral  SnhmndPl^ 
Material s 

An  empty  stainless  steel  HPLC  column,   50.0  cm  long  with 
a  2.54  cm  outside  diameter,  and  2.25  cm  inside  diameter  was 
used  for  column  experiments.  Columns  of  this  size  are 
typically  used  as  preparative  HPLC  columns.  The  column 
packing  was  a  1:250   (w:w)  mixture  of  AG  MP-50  Analytical 
grade  Macroporous  Cation  Resin  and  pure  quartz  sand  derived 
from  a  soil  sample  removed  from  the  E2  horizon,   at  a  depth  45 
to  70  cm  in  a  soil  classified  as  a  St.  Lucy  fine  sand 
(Location  of  pit  NW  1/4  of  SE  1/4  Sec.  T32S  R20E) .  The  soil 
sample  was  sieved  and  any  organic  material  was  oxidized  by 
repeated  washings  with  30%   (volume: volume)  H2O2.  Mariotte 
bottles  were  made  from  1  dm^  Wheaton  "4  00"  brand  borosilicate 
glass  serum  bottle.  A  Rainin  Rabbit  Peristaltic  pump  (Rainin 
Instrument  Co.,   Inc.;Mack  RD;  Woburn,  MA  01801)  mounted  with 
Fisherbrand  PVC  Manifold  Tubing  with  a  1.85  mm  inside 
diameter  was  used  to  drive  the  solutions  through  the  packed 
column.  Aside  from  the  manifold  tubing,   all  tubing  was  Tygon™ 
Clear  laboratory  tubing   (3.2  mm  inside  diameter,   6.4  mm 
outside  diameter) .  Choice  of  column  influent  was  made  by 
using  a  Pharmacia  solvent  resistant  valve  SRV- 3 (supplied  by 
Pharmacia  Inc.;   800  Centennial  AV;  Piscataway,  NJ  08854).  The 
flow  rate  was  always  0.0167  g  s~l   (1  g  min"!) .  This 
corresponds  to  a  nominal  linear  flow  rate  of  4.20x10-5  m  s"! 
and  an  interstitial  velocity  of  1.31x10-^  m  s"!. 
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Experimental  Procedure 
Determination  nf  colnmn  CKC 

The  column  was  saturated  with  0.01  mol   (NaC104)  dm~3 
solution.  The  mass  of  the  column  was  determined.  The  influent 
solution  was  changed  to  a  0.01  mol   (1/2  Ca(C104)2)  dm-3 
solution.  Fractions  of  the  column  effluent  were  collected. 
Sodium  in  these  factions  was  determined  by  flame  emission 
spectrometry.  The  column  CEC   (mol   (Nal+) )  was  taken  to  be  the 
total  Na  eluted  (mol  Nal+)   less  the  Na  in  solution  in  column 
before  displacement  by  Ca2+  (mol  Na^"*")  . 

Determination  of  mass  of  solids 

The  empty  column  including  the  tube,  frits,  ferrules, 
end  fitting  bodies,  nuts  and  plastic  male  column  end  caps, 
was  weighed.  The  packed  column  was  oven-dried  until  the 
measured  mass  of  the  packed  column  reached  a  constant  (to  a 
tenth  of  a  gram)  value.  The  mass  of  solids  was  taken  to  be 
the  mass  of  oven-dry  packed  column  less  those  of  the  empty 
column  and  the  end  fitting  bodies 

Determination  of  porosity 

As  noted  above,  the  majority  of  the  column  packing  was 
quartz  sand.  Wilding  et  al.(1977)   cited  a  value  of  2.65  kg 
dm-3  for  the  density  of  quartz.  The  porosity  (Scheidegger, 
1972)  was  therefore 


48 


porosity  =  (total  volume  -  (particle 

weight/particle  density) ) /total  volume 

Determination  of  vnlnmPtrin  water  contpnl- 

The  mass  of  the  water-filled  packed  column  was 
determined  after  every  run.  The  volumetric  water  content  was 
calculated  by 

Volumetric  water  content  =  Total  volume  - 

(  (Total  mass  -  mass  column  - 
mass  solids) /density  of  water) 

■Column  ExperimAntf? 

All  column  experiments  were  conducted  in  a  laboratory 
with  ambient  temperature  at  24 . 5  ±  0 .  6  C° .  A  schematic  of  th 
arrangement  of  the  equipment  used  in  the  column  experiments 
is  presented  in  Figure  3. 

First  case,  R^3n  A.  All  solutions  passing  through  the 
column  were  saturated  with  He  gas  to  reduce  the  formation  of 
air  bubbles  in  the  column  by  cavitation.  At  least  15  pore 
volumes  of  0.01  mol(l/2  Ca (0104)2)   dm-3  solution  were  passed 
through  the  column  with  a  flow  rate  of  1  g  min"!-   The  pump 
was  stopped  momentarily  so  that  the  influent  stream  could  be 
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changed  to  a  solution  containing  0.005  mol   (1/2  Mg(C104)2) 
dm-3  0.005  mol (NaC104 ) dm-3  solution.  After  61  minutes,  the 
influent  flow  was  reverted  to  0.01  mol  (1/2  Ca  (0104)2)  dm-3 
solution.  Fifty-five  fractions  of  the  effluent  of 
approximately  10  g  were  collected.   The  solution 
concentrations  of  Ca,  Mg,   and  Na  were  determined  for  each 
effluent  sample  by  the  methods  described  in  the  section  of 
this  chapter  describing  the  batch  exchange  determinations. 

Second  cgsSr  Run  P.  All  solutions  passing  through  the 
column  were  saturated  with  He  gas  to  reduce  the  formation  of 
bubbles  in  the  column  by  cavitation.  At  least  15  pore  volumes 
of  0.01  mol(NaC104)   dm-3  solution  were  passed  through  the 
column  at  a  flow  rate  of  1  g  min-L  The  pump  was  stopped 
momentarily  so  that  the  influent  stream  could  be  change  to  a 
solution  containing  0.005  mol (1/2  Ca (CIO4) 2) dm-3  q.OOS  mol (1/2 
Mg(C104)2)   dm-3  solution.  After  61  minutes,   the  influent  flow 
was  reverted  to  0.01  mol(l/2  Ca(C104)2)   dm-3  solution. 
Fifty-nine  fractions  of  the  effluent  of  approximately  10  g 
were  collected.   Solution  concentrations  of  Ca,  Mg,   and  Na 
were  determined  by  the  methods  described  in  the  section  of 
this  chapter  describing  the  batch  exchange  determinations. 


COMPUTER  MODEL  DEVELOPMENT 


Two  computer  programs  were  used  in  the  evaluation  of 
three  exchanger  phase  reference  functions  for  inclusion  as 
chemical  submodels  of  cation  transport  models.  The  first 
computer  program,  named  BLCKBX  (Black  box) ,  calculated  the 
binary  and  ternary  isotherms  dictated  by  the  thermodynamic 
exchange  constants,   solution  phase  and  exchanger  phase 
activity  coefficients.  The  exchange  isotherms  calculated  by 
Program  BLCKBX  were  used,   along  with  appropriate  initial  and 
boundary  conditions,   in  the  second  program,  TABMODEL  (Tabular 
model),  to  calculate  cation  transport.  Two  separate  programs, 
rather  than  one  integrated  program,  were  developed  because  it 
was  found  to  be  more  efficient  computationally  to  calculate  a 
complete  exchange  isotherm  before  attempting  to  solve  the 
solute  transport  problem.  The  development  and  uses  of  these 
programs  will  be  discussed  below. 


The  Comonter  Program  RT.PKRX 


Program  BLCKBX  is  a  FORTRAN  computer  program  developed 
to  calculate  instantaneous  partitioning  of  three  cations  of 
known  amounts  between  the  exchanger  and  solution  phases.  A 
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listing  of  this  program  is  presented  as  an  appendix  to  this 
dissertation.  The  program  was  used  to  calculate  equilibrium 
exchange  isotherms  to  be  used  to  represent  chemical  submodels 
in  convective-dispersive   (advective-dispersive)  transport 
models  of  systems  involving  steady  state  fluid  flow  through 
porous  media. 

The  program  BLCKBX  calculated  ternary  exchange 
equilibria  constrained  by  charge  and  mass  balance  in  the 
exchanger  and  solution  phases.  Briefly,   for  each  iteration, 
one  of  the  exchange  selectivity  coefficients  was  used  to 
estimate  the  distribution  coefficient  of  one  of  the  cations. 
The  distribution  coefficient    (di  )   is  the  ratio  of  a  cation  in 
the  exchanger  phase  to  its  concentration  in  the  solution 
phase  or 


e 
m. 


~     3  [40] 


where 

s  _ 

-  cation  i  solution  phase  concentration  (mol  dm-3) , 

and 

e  _ 

m.  =  cation  i  exchanger  phase  molality   (mol  kg-l)  . 

Since  the  total  amount  was  specified  as  input,  the 
solution  phase  concentration  and  exchanger  phase  molality  of 
the  cation  in  the  two  phases  could  be  calculated  directly. 
The  solution  phase  concentration  and  exchanger  phase  molality 
of  the  other  cation  in  the  exchange  selectivity  coefficient 
were  calculated  by  difference.  These  estimates  were  then  used 
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to  calculate  the  distribution  coefficient  in  the  next 
iteration.  The  program  allowed  the  selection  of  one  of  three 
exchanger  phase  reference  functions:  the  mole  fraction  used 
in  the  Vanselow  (1932)   formulation,  the  equivalent  fraction 
found  in  the  work  of  Gaines  and  Thomas   (1953),  and  the 
statistical  thermodynamic  suggested  by  Krishnamoorthy  et  al. 
(1948)  . 


Input  to  Prnaram  BLCKRX 

An  example  of  an  input  file  for  program  input  is 
included  as  an  appendix  to  this  dissertation.  Aside  from 
variables  which  specify  the  iterative  tolerances  and 
iterative  limits  of  the  computer  program,   seven  variables  or 
groups  of  variables  were  specified  for  each  implementation  of 
the  program: 

1.  Porous  medium  bulk  density   (BULKD)   in  kilograms  per 
cubic  decimeter; 

2.  Volumetric  water  content   (THETA)   in  cubic  decimeters 
per  cubic  decimeter; 

3.  Cation  exchange  capacity  (CECFIX)  in  moles   (+)  per 
kilogram; 

4.  The  sum  charge  born  by  all  cations  in  solution 
(CTOTAL)   in  moles (+)  per  cubic  decimeter; 

5.  The  three  thermodynamic  exchange  constants  (RKIJO, 
RKJKO,and  RKKIO;   corresponding  to  i^exi/j,  i<exjyjt,  and  i?exit/i, 
respectively)   which  were  dimensionless ; 

6.  The  valences   (ZI,   ZJ,   and  ZK)   of  each  cation; 
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7.  The  total  moles  of  each  cation  in  the 
exchanger-solution  system  (CCI, CCJ, CCK)   in  moles  per  cubic 
decimeter . 

inclusion  of  variable  SelerMvities  anH  ActivitY  Hnef  f  i  ri  Pnt 

In  general,  exchange  select ivities  are  not  fixed,  but 
vary  with  changes  in  the  composition  of  the  exchanger  phase. 
Estimates  of  these  changes  could  be  included  in  the 
calculation  of  exchange  equilibria  by  several  means.  Three 
subroutines  could  be  used  to  introduce  estimates  of  exchanger 
phase  and  solution  phase  activity  coefficients  or  exchange 
select ivities  which  vary  with  changes  in  the  exchanger  phase 
or  solution  phase  composition.  The  exchanger  phase  and 
solution  phase  coefficients  were  called  Fl,  fj,  or  FK  and 
GAMMAI,   GAMMAJ,    or  GAMMAK,   respectively.   Every  time  the 
defining  relations  were  used  in  a  calculation,  the  values  of 
the  exchanger  phase  or  solution  phase  activity  coefficients 
could  be  adjusted  by  calling  SUBROUTINE  F,   SUBROUTINE  GAMMA, 
and  SUBROUTINE  KCALC.   SUBROUTINE  KCALC  calculated  the 
exchange^^lej::tiyity  c  if  the  exchanger  phase 

activity  coefficients  were  not  estimated,  SUBROUTINE  KCALC 
could  be  used  to  insert  exchange  selectivity  coefficient 
values  estimated  by  least-squares  regression  of  the 
relationship  between  the  ion  exchange  selectivity 
coefficients  to  the  exchanger-solution  system  composition. 
If  one  had  a  regression  equation  relating  the  natural 
logarithm  of  the  Vanselow  selectivity  coefficient.   In  k^,,,/j, 
and  the  exchanger  phase  cation  i  equivalent  fraction  and  a 
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second  regression  equation  relating  In  A:(v)j/;c  to  exchanger 
phase  cation  k  equivalent  fraction,   the  third  ion  exchange 
selectivity  coefficient,   A:(v)Vi  '   could  be  calculated  by  the 
triangle  rule,   one  form  of  which  states  that 

Zklnk^^)i/j  +  zilnk(y)  j/i^  +  ZjlnA:(v)  jt/i  =  0  [41] 

The  three  estimates  of  the  ion  exchange  selectivity 
coefficients  could  be  calculated  within  SUBROUTINE  KCALC, 
each  iteration  of  the  main  program  being  adjusted  as  the 
estimated  exchanger  phase  composition  was  changed.  While  all 
implementations  of  variable  exchange  select ivities  described 
in  this  study  used  a  variable  adjusted  Vanselow  selectivity 
coefficient,  that  is  the  value  of  the  adjusted  Vanselow 
selectivity  coefficient  varied  as  a  function  of  the  exchanger 
phase  composition;   in  principle,   the  same  results  could  be 
achieved  using  other  variable  select ivities .  The  program 
presented  the  solution  phase  concentrations  and  exchanger 
phase    molalities   (in  moles  per  decimeter  and  moles  per 
kilogram,   respectively)   as  output. 
The   sSteps   in  the  Program  BLCKRX  Algorithm 

The  steps  in  the  algorithm  are  given  below. 
Step  1.  Given  the  ratios  of  cation  exchange  capacity 
(CEC)   and  charge  born  by  all  cations  in  solution,   the  program 
gave  first  estimates  of  the  exchanger  phase  and  solution 
phase  concentrations. 
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CEC 

di(est)=  — -  [42] 


3 


where 


s 


Cj-^  =  sum  of  the  charge  born  by  all  cations  in  solution 

(mol(  +  )  dm-3) 

Step  2.  Based  on  the  estimated  distribution 
coefficients,   the  program  calculated  initial  values  of  the 
exchanger  phase  and  solution  phase  concentrations  of  all 
cations . 


s  _  Nj 


_e      ^  Ni 

9 


~     0  [44] 


where  Ni  =  total  concentration  of  cation  solute  i  in  both  the 
exchanger  and  solute  phases    (mol  dm~3) 


Ni  =  Pb        +  0  cl 


[45; 


Step  3.  Using  the  initial  values  of  the  exchanger  phase 
and  solution  phase  concentrations,  the  program  recalculated 
the  distribution  coefficient  for  one  of  the  ions 


^    -  — -  [4  6- 

liCi 


57 


Step  4.  Based  on  the  estimated  distribution 
coefficients,  the  program  calculated  initial  values  of  the 
exchanger  phase  and  solution  phase  concentrations  of  all 
cations 


3   Ni 

c.  =   ^ 


Ph  di  +0 


m.  = 


e 


Step  5.  Using  these  exchanger  phase  and  solution  phase 
concentrations,   the  program  calculated  by  difference  the 
exchanger  phase  and  solution  phase  concentrations  of  the 
other  cation  in  the  defining  relation. 


CEC  -  /,Zkml-Zim^. 
k9ii ;  j 


e 

777  .  = 


Zj  [47 


s 

c  .  = 


k?ti ;  j 


^  Zj  [48] 


Step  6.   The  program  then  checked  to  determine  if  both 
s 

ij    and  Cj    were  greater  than  zero. 

If  the  two  were  greater  than  zero.   Steps  3-6  were 
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repeated  until  the  estimated  exchanger  phase  composition 
converged  to  an  acceptable  absolute  difference  between 
iterations.    (Please  note  that  during  the  solution  of  the 
exchanger  phase  composition  using  i<exi/j  the  concentrations  of 
all  cations  other  than  cation  i  and  cation  j  were  fixed;  that 
is,   the  sums 

e 

^k^k  I  and 
k9ii ;  j 

n 

k;ti ;  j 


in  equations   [47]  and   [48]   were  fixed.  The  other  cations  were 
allowed  to  vary   (and  thus  have  their  exchanger  phase  and 
solution  phase  concentrations  estimated)   by  following  the 
same  Steps  1-6  but  using  with  different  thermodynamic 
exchange  constants    (e  .  g  .  i^exj/;t  or  K^^^l  i  )•) 

If  one  of  the  two  was  not  greater  than  zero,  the 
program  assigned        ,   r^.  ,   c\  ,  and    rn-     their  previous 
values  and  used  the  same  defining  relationship  to  estimate 
dj   (instead  of  di)   and  the  Steps  4-6  were  carried  out  with 
roles  of  cation  i  and  cation  j  reversed.   If,  at  Step  6,  c\ 
and  m\    were  not  greater  than  zero  the  program  went  to  Step  7 
and  selected  another  defining  relation   (e.g.  i^exj/;c)  to 
estimate  iteratively  the  exchanger  phase  composition. 
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For  an  n-ary  system  (i.e.  one  with  n  cations)   there  are 
2*n  possible  iterative  loops;  only  n-1  iterative  loops  need 
converge  to  estimate  the  exchanger  phase  composition. 

Step  7.  Once  the  estimates  for  the  defining  relation 
involving  /Cexi/j    have  been  resolved,   the  program  performed 
Steps  1-6  using  another  ion  exchange  constant.  As  before,  the 
program  checked  for  convergence  and  non-negative  estimated 
exchanger  phase  and  solution  phase  concentrations. 

Step  8.  The  program  repeated  Step  7  until  n-1  defining 
relations  had  been  resolved. 

Step  9.  The  program  used  the  resulting  estimates  of 
solution  phase  and  exchanger  phase  compositions  as  the  input 
for  a  repetition  of  Steps  3-8  until  the  estimated  exchanger 
phase  composition  converged  to  an  acceptable  absolute 
difference  between  iterations. 

For  reproducible  estimates,  meaning  predicted  isotherms 
which  are  identical  regardless  of  the  order  in  which  the 
defining  relations  were  chosen  in  the  algorithm,  the 
tolerance  used  to  judge  convergence  of  exchanger  phase 
composition  had  to  be  on  the  order  of  lO'S.  This  caused 
millions  of  iterations  to  be  performed  so  that  the  resolution 
of  a  ternary  isotherm  consisting  of  2500  distinct  fractions 
of  cationic  amounts  could  take  1.5  hours  of  central 
processing  unit    (CPU)   time  on  a  Digital  Equipment  Corporation 
VAX  11/750  minicomputer  having  a  floating  point  coprocessor. 
No  attempt  to  optimize  execution  time  while  preserving 
reproducibility  was  attempted. 
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The  accuracy  of  BLCKBX  as  a  method  of  calculating 
ternary,  heterovalent  exchange  equilibria  was  checked  by 
three  tests.  First,  the  output  of  the  program  maintained 
overall  cation  mass  balance  between  the  exchanger  and 
solution  phases,   and  maintained  charge  balance  in  both  the 
solution  and  exchanger  phases.  Second,  the  calculated  cation 
exchanger  phase  molalities  and  cation  solution  phase 
concentrations  consistently  reproduced  each  of  the 
appropriate  exchange  selectivity  coefficients.  Third,  program 
BLCKBX  generated  the  same  exchange  isotherm  irrespective  of 
the  order  by  which  the  defining  relations  were  solved  within 
the  program. 


The  Computer  Program  TABMODEL 


The  second  computer  program,   TABMODEL,  used  the 
isotherms  generated  by  BLCKBX  to  solve  the  simultaneous 
one-dimensional  transport  of  three  cationic  solutes  through 
porous  media.  The  program  was  developed  by  S.A.  Bloom  and 
R.S.  Mansell  of  the  Soil  Science  Department,  University  of 
Florida.  The  program  TABMODEL    used  a  tabular  look-up  method 
of  solving  the  one-dimensional  convective  dispersive 
equation  for  multicomponent  cation  transport  through  porous 
media.  The  convective-dispersive  (advective-dispersive) 
equation 
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was  approximated  numerically  according  to  the 
Crank-Nicholson  convention  by 


c/,  (z,n  +1)   -  c^,(z,n) 


At  D 
2Az  2 


[c/,{z-l,r!  +1)   -  2c/,(z,/i  +1)   +     c/,(z  +  l,n 


+1) 


+  c/,  (z-l,n)  -  2c/,  (z,n)   -  c/,  (z  +  ] 


At  V 
4Az 


[c/.(z  +  l,n  +  1)   -  c/,(z-l,n  +1)   +  c/,(z  +  l,n) 


-    c/,  (Z-1,  71)  ] 

"  ~^  +1)   -m^,(z,n))  [49] 


^     ,    e  .  .  .  e 


By  collecting  terms,  this  could  be  represented  by 


At  D       At  V 
?Az  2      4Az  • 
g       At  v-i  s  r     At  D 


+  c^,(z+l,n  +1    I   +   J  = 

2Az  2  4Az 


3  ,      ^  r  At    D         At    VT  »  r  I 

CiAz-l,n)  l—^  +  ——J    +  c/,(z,n  )     [  1-   , 

2Az  4Az  ^  Az  2 

s  r  At  D       At  VI  Pk 

.  c,.,z.i,n,  [—  -—]-f-  .;,,z,.  .1, 


[50] 
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As  in  the  example  in  the  literature  review,  this 
formulation  was  in  turn  converted  to  matrix  form  and  solved 
by  matrix  multiplication. 

This  algorithm  differed  from  other  approaches  in  that 
the  solution  of  the  solution  phase  concentration  depended  on 
the  final  equilibrium  concentration.  There  were  four  basic 
steps  in  the  solution  of  the  algorithm: 

1.  Calculate  estimate  of  c/,  (z,n  +1)  based,   in  first 
iteration  in  a  loop  m^,(z,n  +1)   =  mf,{z,n)   is  accepted; 

in  later  iterations  the  estimate  of  m^,{z,n  +1)  was  based  on 
the  estimate  resulting  from  the  previous  iteration; 

2.  Using  the  solution  and  exchanger  phase 
concentrations  the  program  calculated  total  cation 
concentration  in  an  element; 

3.  Using  an  isotherm  calculated  in  advance  by  program 
BLCKBX,  the  program  calculated  equilibrium  distributions  of 
the  three  cations  in  the  exchanger  and  solution  phases; 

4.  The  program  then  used  the  resulting  estimate  of 
exchanger  phase  composition  to  present  a  new  estimate  of 
m^,(z,n  +1)  . 

The  algorithm  used  the  new  estimate  of  exchanger  phase 
composition  and  performed  steps  1-4  until  convergence.  The 
user  of  the  computer  program  was  free  to  choose  any  method 
for  calculating  the  equilibria  which  gave  the  estimates  of 

m^Az,n  +1).  This  freedom  allowed  the  evaluation  of  both 
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variable  selectivity  and  several  exchanger  phase  reference 
functions  for  the  prediction  of  cation  transport  in  ternary 
systems . 

As  will  be  shown  later,  the  combination  of  these  two 
computer  programs  was  able  to  duplicate  accurately 
experimental  column  data  and  the  simulations  based  on  an 
algorithm  suggested  by  Valocchi  et  al.  (1981b). 


RESULTS  AND  DISCUSSION 
Batch  Expprimpnf.q 

Measured  supernatant  solution  concentrations  and 
exchangeable  cation  contents  from  all  replicates  of  the 
exchange  experiments  with  the  cation  resin  are  presented  in 
Tables  3  through  9.  Mean  values  of  supernatant  solution 
concentrations  and  exchangeable  cation  mole  fractions  for 
each  equilibration  composition  are  presented  in  Tables  10  and 
11.  These  mean  data  were  used  in  all  subsequent  analyses  of 
the  batch  exchange  data.  The  data  reported  by  Elprince  et  al. 
(1980),  who  measured  NH4I+-  Ba2+-La3+  exchange  equilibria  on 
montmorillonite  shaken  with  salt  solutions  having  constant 
0.1  mol  dm-3  ci  concentrations,   were  used  to  compare  ternary 
exchange  on  natural  exchangers  with  ternary  exchange  on 
synthetic  exchangers, 

VanSfilow  COeffirjent  on  exrhanapr  phase  egniv^l^nt 
fraction  regression 

The  adjusted  Vanselow  selectivity  coefficients  were 
calculated  for  binary  and  ternary  exchange  experiments  with 
the     cation  resin  and  those  with  montmorillonite  reported  by 
Elprince  et  al.    (1980).     For  both  groups  of  binary  exchange 
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Table  7.  Experimental  exchange  data  for  Set  Five 


Supernatant  Exchangable 
cation  concentrar  1nn  cation  rnnfAnt- 


Equilibration     Replication      Ca  Mg  Ca 


Mg 


A 

1 

A 

2 

A 

3 

B 

1 

B 

2 

B 

3 

C 

1 

C 

2 

C 

3 

D 

1 

D 

2 

D 

3 

E 

1 

E 

2 

TP 

hi 

3 

F 

1 

F 

2 

F 

3 

G 

1 

G 

2 

G 

3 

H 

1 

H 

2 

H 

3 

I 

1 

I 

2 

I 

3 

J 

1 

J 

2 

J 

3 

K 

1 

K 

2 

--iTunol 

dm 

_-3 
o  

c 
O 

A  n 

.  028 

0 

.000 

c 

5 

A  T 

.087 

0 

.  000 

c 
D 

.087 

0 

.  000 

0 

.  162 

4 

.440 

0 

.173 

4 

.  407 

0 

.  170 

4 

.351 

0 

.  964 

4 

.  144 

1 

.  006 

4 

.110 

0 

.  996 

4 

.092 

1 

.471 

3 

.  695 

1 

.474 

3 

.  644 

1 

.  452 

3 

.  585 

1 

.  958 

3 

.  187 

1 

A  O  T 

.  927 

3 

.  160 

006 

3 

201 

U 

2 

683 

U 

DUO 

2 

729 

U 

000 

2 

706 

o  "3 

yjo 

2 

132 

2. 

867 

2  . 

132 

2. 

880 

2  . 

081 

3. 

417 

1. 

637 

3. 

369 

1 . 

643 

3. 

327 

1 . 

643 

3. 

840 

1 . 

189 

3. 

933 

1. 

130 

3. 

865 

1 . 

130 

4. 

332 

0. 

565 

4  . 

216 

0. 

571 

4  . 

139 

0. 

568 

4  . 

759 

ND 

4  . 

811 

ND 

mmol 


0 

.344 

0 

.000 

0 

.293 

0 

.000 

0 

.303 

0 

.000 

0 

.116 

0 

.216 

0 

.120 

0 

.219 

0 

.  119 

0 

.218 

0 

.249 

0 

.180 

0 

.232 

0 

.168 

0 

.229 

0 

.172 

0 

.256 

0 

.111 

0 

.262 

0 

.124 

ND 

ND 

0 

.304 

0 

.094 

0 

.296 

0 

.086 

0 

.289 

0 

.078 

0 

361 

0 

.052 

0 

274 

0 

045 

0 

293 

0 

044 

0 

292 

0 

043 

0 

340 

0 

049 

0. 

353 

0 

050 

0. 

291 

0. 

027 

0. 

362 

0. 

036 

0. 

329 

0. 

033 

0. 

353 

0. 

021 

0. 

335 

0. 

021 

0. 

321 

0. 

022 

0. 

309 

ND 

0. 

270 

0. 

008 

0. 

364 

ND 

ND 

ND 

ND 

ND 
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Table  8.  Experimental  exchange  data  for  Set  Six. 


Supernatant  Exchangable 
cation  concent rrit ion  cation  r-nni-Ani- 


Equilibration    Replication      Ca  Na  Ca 


Na 


A 

1 

A 

2 

A 

3 

B 

1 

B 

2 

B 

3 

C 

1 

C 

2 

c 

3 

D 

1 

D 

2 

D 

3 

E 

1 

E 

2 

E 

3 

F 

1 

F 

2 

F 

3 

G 

1 

G 

2 

G 

3 

H 

1 

H 

2 

H 

3 

I 

1 

I 

2 

I 

3 

J 

1 

J 

2 

J 

3 

K 

1 

K 

2 

 mmol  dm~3  


0 

.021 

10.08 

0 

.016 

10.24 

0 

.015 

10.03 

0 

.003 

9.889 

0 

.000 

10.12 

0 

.016 

9.936 

0 

.065 

9,  936 

0 

.  041 

9.843 

0 

.047 

9.982 

0 

.133 

9.  936 

0 

.154 

9.750 

0 

.139 

9.  982 

0 

.285 

10.01 

0 

266 

9.811 

0 

304 

9.  962 

1 

148 

9.261 

0 

636 

9.311 

0. 

650 

9.4  60 

0. 

887 

8.764 

0. 

889 

8.764 

0. 

901 

8.764 

1 . 

949 

6.932 

1 . 

948 

6.764 

1. 

943 

6.731 

3. 

215 

3.748 

3. 

270 

3.816 

3. 

300 

3.748 

4  . 

546 

3.475 

4  . 

517 

3.398 

4. 

595 

3.411 

4  . 

947 

0.000 

4  . 

828 

0.000 

mmol 


0.000 

0 

.791 

0.000 

0 

.732 

0.000 

0 

.837 

0  .114 

0 

.407 

0.137 

0 

.506 

0.129 

0 

.485 

0.234 

0 

.376 

0.212 

0 

.327 

0.237 

0 

.370 

0.321 

0 

.234 

0.278 

0 

.237 

ND 

ND 

0.345 

0 

.  138 

0.339 

0 

.  153 

0.337 

0 

108 

0.351 

0 

064 

0.362 

0 

070 

0.335 

0 

053 

0.337 

0. 

034 

0.358 

0. 

035 

0.338 

0. 

052 

0.411 

0. 

026 

0.384 

0. 

031 

0.301 

0. 

015 

0.349 

0. 

003 

0.334 

0. 

009 

0.344 

0. 

012 

0.207 

ND 

0.392 

0. 

000 

0.379 

ND 

0.362 

0. 

003 

0.343 

0. 

004 
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Table  9.  Experimental  exchange  data  for  Set  Seven. 


Supernatant  Exchangable 
cation  roncenl-rrir i nn  cation  mnt-^n^- 


Equilibration     Replication      Mg  Na 


Mg  Na 


A 

1 

■3 

R 

1 

1 

o 

R 

o 

r 
\^ 

1 

1 

r 

o 

J 

n 
u 

1 

n 

Z 

n 

o 
J 

F 

Hi 

1 

Hi 

id 

Hi 

3 

F 

1 

F 

2 

F 

3 

G 

1 

G 

2 

G 

3 

H 

1 

H 

2 

H 

3 

I 

1 

I 

2 

I 

3 

J 

1 

J 

2 

J 

3 

K 

1 

K 

2 

 mmol  dm  ^  


U  .  UUU 

9 .  991 

U  .  UUU 

10.04 

U  .  UUU 

9 .  941 

u .  u  i  y 

10.09 

U  .  U  1  D 

10.14 

V  »\J  1  / 

10.14 

u .  Uo  J 

9 .  941 

A    A  T  £r 
U  .  U  /  D 

9.741 

A  ATI 

U  .  U71 

9.642 

A     T  1  A 

U  .  IIU 

9.891 

A      1  A  O 

U  .  103 

9.691 

U  .  lU  J 

9.791 

A     A  yl  A 
U  .  U4  U 

9.941 

A     A  /I  O 

(J  .  U43 

9 .  941 

A     A  O  1 

U  .  (J37 

9.841 

0.261 

9   4  4  R 

0.252 

9.345 

0.274 

9.396 

0.483 

9.139 

0.396 

9. 191 

0.402 

9. 139 

1 .726 

6.804 

1.722 

6.597 

1.737 

6.597 

3.191 

3.709 

3.137 

3  .815 

3.164 

3.744 

4  .355 

1 . 128 

4  .400 

1 .115 

4  .355 

1 .117 

4  .  917 

0.000 

4  .850 

0.000 

mmol 


0 

.  000 

0 

.812 

0 

.  000 

0 

.866 

0 

.  000 

0 

.885 

0 

.042 

0 

.736 

0 

.043 

0 

.720 

0 

.  044 

0 

.706 

0 

.  048 

0 

.632 

0 

.051 

0 

.  994 

0 

.  048 

0 

.629 

0 

.  152 

0 

.432 

0 

.  162 

0 

.450 

0 

.165 

0 

.506 

0 

.  120 

0 

.542 

0 

.135 

0 

.532 

0 

.  129 

0 

.572 

0 

.266 

0 

.282 

0 

.264 

0 

283 

0 

243 

0 

283 

0 

304 

0 

197 

0 

301 

0 

227 

0. 

303 

0. 

234 

0. 

351 

0. 

084 

0. 

363 

0. 

086 

0. 

334 

0. 

075 

0. 

389 

0. 

041 

0. 

389 

0. 

042 

0. 

381 

0. 

043 

0. 

342 

0. 

Oil 

0. 

351 

0. 

012 

0. 

348 

0. 

Oil 

0. 

358 

0. 

004 

0. 

343 

0. 

005 
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Table  12.  Binary  exchange  equilbria  with  montmorillonite 
reported  by  Elprince  et  al .  (1980)  . 


Cation  .solution  conrgntcvati on    Exch;^na;.hie  cation  nn^f^^^ 


La 


Ba 


NH4 


La 


Ba 


NH4 


0.00350 
0.00524 
0.00669 
0.00783 
0.00863 
0.00175 
0.01575 
0.02075 
0  .02671 

0.02618 
0.02074 
0.01599 
0.01191 
0.00859 
0.00764 
0.00641 
0.00503 
0.00307 


-mol  dm~3- 

0.04455 
0.04170 
0.03945 
0.03776 
0.03665 
0.03183 
0.02573 
0.01832 
0.00959 


0.04103 
0.03391 
0.02832 
0.02425 
0.02119 
0.01868 
0.01586 
0.01249 
0.00774 


0.01914 
0.03571 
0.05071 
0.06386 
0.07314 
0.07600 
0.07940 
0.08343 
0.08814 

0.0814 

0.03271 

0.04443 

0.05314 

0.05886 

0.06386 

0.06900 

0.07543 

0.08457 


0. 1560 
0.3987 
0. 6543 
0. 9253 
1 .2140 
1 .0560 
0.8550 
0. 6073 
0.3190 

0.296 
0.580 
0.826 
1.244 
1.446 
1.108 
0.  657 
0.398 
0.169 


mmol- 


0.2595 
0.4015 
0.5255 
0  .  6345 
0.7340 
0 .4400 
0.2280 
0.0895 
0.0220 


0.153 

0.327 

0.565 

1  .085 

2  .221 

1.536 

1 .094 

0.871 

0.583 

9.05 

1.9 

16.26 

9.38 

21 .  96 

16.08 

26.12 

27.77 

29.24 

42.3 

21.  6 

37.1 

14.3 

31.7 

7.59 

25.0 

2.29 

15  .  68 
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Table  13.  Ternary  exchange  equilbria  with  montmorillonite 
reported  by  Elprince  et  al .  (1980)  . 


Cation 

solution 

_QOncentratiori 

Exchanaable  ration 

contpnt 

La 

Ba 

NH4 

La  Ba 

NH4 

0.00138 
0.00181 
0.00200 
0.00216 
0.00302 
0.00442 
0.00738 
0.01142 


0.00310 
0.00298 
0.00640 
0.00770 
0.00075 
0.00299 
0.01233 
0.02033 


-mol  dm~3- 

0.04321 
0.03866 
0.03499 
0.03202 
0.02674 
0.02283 
0.01708 
0.00940 

0.00884 
0.01589 
0.02367 
0.03038 
0.02219 
0.02175 
0.01775 
0.00955 


mmol- 


0 

.0096 

0 

.1933 

0 

.  655 

0 

.120 

0 

.0171 

0 

.4810 

1 

.117 

0 

.449 

0 

.0237 

0 

.7923 

1 

.511 

0 

.857 

0 

.0294 

1 

.  1063 

1 

.863 

1 

.352 

0 

.0372 

1 

.3500 

1 

.407 

1 

.547 

0 

.0404 

1 

.2083 

0 

.755 

1 

.149 

0 

0440 

0 

9120 

0 

.307 

0 

.716 

0 

0473 

0 

5107 

0 

063 

0 

.319 

0 

0732 

0 

3523 

0 

1205 

0 

712 

0 

0591 

0 

6950 

0 

4295 

1 

185 

0. 

0331 

1 . 

3471 

0. 

6760 

0, 

899 

0. 

0158 

1 . 

8673 

1. 

0400 

0. 

747 

0. 

0540 

0. 

5873 

1. 

8815 

2  . 

824 

0. 

0475 

1 . 

0237 

0. 

8720 

1. 

424 

0. 

0271 

1 . 

0753 

0. 

2400 

0. 

422 

0. 

0194 

0. 

6060 

0. 

0485 

0. 

127 
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data,   least  squares  regressions  were  performed  to  fit  the 
data  to  the  following  equation 

where  Xj-    is  the  exchanger  phase  equivalent  fraction  of 
cation  "j".     The  parameter  estimates  and  coefficients  of 
determination  are  presented  in  Tables  14  and  15.   In  one  of 
the  binary  exchange  experiments  with  montmorillonite,  the 
parameter  p2  could  be  dropped  from  the  regression  of 
-^(vjLa/NH^        exchanger  phase  NH4  equivalent  fraction 

without  decreasing  the  coefficient  of  determination.  The 
thermodynamic  exchange  constants  were  calculated  using  the 
regression  equations  by  evaluating  the  integral 


1 

In  Ke^i/j  = 

0 


that  is 


The  resulting  estimated  exchange  constants  were  tested  by  the 
triangle  rule 

Zk  Ini^exi/j  +  Zi  In/Ce^j/jt  +  Zj  lnfCex;c/i  =  0. 
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The  parameters  for  each  of  the  three  quadratic  regression 
parameters  were  adjusted  so  that  the  triangle  rule  was  obeyed 
exactly.  These  adjusted  parameter  estimates  are  presented  in 
Tables  16  and  17.  The  thermodynamic  exchange  constants  and 
standard  free  energies  of  exchange  for  the  cation  resin  and 
montmorillonite  calculated  using  the  adjusted  regression 
parameters  are  presented  in  Table  18.  The  standard  free 
energies  of  exchange  measured  in  this  study  had  higher 
absolute  value  than  those  reported  by  Boyd  (1970) .  The 
difference  was  greatest  for  the  Mg/Ca  exchange  reaction. 
These  differences  may  be  due  to  differences  in  electrolyte 
solution  concentrations,  methods  of  data  analysis,  or 
specific  effects  due  to  the  common  anion. 

The  regression  equations  and  natural  logarithms  of  the 

adjusted  Vanselow  selectivity  coefficients  determined  in  both 

the  binary  and  ternary  exchange  experiments  with  the  cation 

resin  and  with  montmorillonite  are  plotted  in  Figures  4 

through  9.  These  plots  show  that,  at  least  qualitatively, 
that  the  relationships  between  In  k^^^^.^.    and  exchanger  phase 
equivalent  fractions  were  similar  in  both  binary  and  ternary 

systems,   though  substantial  deviations  between  the  binary  and 

ternary  data  were  observed.   The  plot  of  In  ;c/,„  against 

{v)Na/Ca  '=1'='-^^'^^ 

exchanger  phase  equivalent  fraction   (Figure  4) determined  in 
ternary  experiments  appeared  to  define  a  different  curve  than 
the  In  values  determined  in  binary  exchange 

experiments.  There  was  so  much  scatter  in  the  plots  of 
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Table  16.  Quadradic  model  parameter  estimates  from  least 
squares  regressions  of  the  natural  logarithm  of  adjusted 
Vanselow  selectivity  coefficients  on  cation  exchanger  phase 
equivalent  fraction  in  binary  systems  with  AG  MP-50 
macroporous  cation  resin  adjusted  so  that  calculated 
thermodynamic  exchange  constants  adhere  to  the  triangle  rule 


 Parameter  est  i  mai-.^,c^ 

Dependant 

Variable 

Po^ 

Pi 

P2 

1^  ■'fivjCa/Mg 

3.447 

-2.665  8 

.  660 

In  A:(v)Mg/Na 

-0. 957 

3.735  -4 

.  672 

1"  ^(v)Na/Ca 

-5.712 

15.784  -12 

.  104 

■•■parameter  estimates  for 

the  model 

1"  ^{v)i/j  =  Po 

+  Pixj .  + 

p2Xj'  2. 

Table  17.  Quadradic  model  parameter  estimates  from  least 
squares  regressions  of  the  natural  logarithm  of  adjusJed 
Vanselow  selectivity  coefficients  on  cation  exchanger  phase 
equivalent  fraction  in  binary  systems  with  montmorinonUe 


 ParamefP]:  (^c-t  t  jr^^r 

Dependant 
Variable 


In 

■''^{v)Ba/La 

-2 

.269 

0.  048 

0  .  910 

In 

-^(v)  La/NH4 

-2 

395 

2  .  076 

In 

^ (v)NH4/Ba 

2 

235 

-3.088 

2.769 

^(v)i/i  =  Po  +  Pixj.    +  ^2Xj.  2. 
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,      ,  a 

^{v)La/NH4  ^'Jamst  exchanger  phase  NH4  equivalent  fraction 
determined  in  both  binary  and  ternary  experiments  that  no 
clear  relationship  between  In  ^c^^^^^a/NH,  exchanger  phase 

NH4  equivalent  fraction  was  apparent,  though  a  linear 

regression  was  used  for  the  purposes  of  this  study.  In  spite 
of  these  exceptions,  the  closeness  of  the  plotted  In 

(v) i/j 

values  against  exchanger  phase  equivalent  fraction 
determined  in  binary  and  ternary  experiments  suggests  that 
relationships  between  In  k^^^^./.    and  exchanger  phase 
composition  determined  in  binary  exchange  experiments  could 
be  valid  frequently  in  predictions  of  ternary  exchange 
equilibria . 

The  results  of  ternary  experiments  were  analyzed 
similarly  to  the  binary  experiment  results.   Two  possible 
regression  models  were  evaluated,   the  simpler  relation  had 
four  parameters 


lnA:(v)i/j  =  Po  +  Pixi    +  P2X2  + 


while  the  more  intricate  model 


[52] 


ln;c(v)i/j  =  Po  +  Pixi    +  P2X2  +  P3X12  +  ^^^^^^  +  p3^^2  f53j 

had  six  parameters.  For  both  equations 

xi  =  exchanger  phase  Na  equivalent  fraction,  and 
X2  =  exchanger  phase  Ca  equivalent  fraction, 
for  the  resin  exchange  data;  and 
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Xi  =  exchanger  phase  NH4  equivalent  fraction,  and 
X2  =  exchanger  phase  Ba  equivalent  fraction 

for  the  montmorillonite  exchange  data  reported  by  Elprince  et 
al.    (1980)  . 

If  equations   [52]  and  [53]  are  to  be  applied 
successfully  to  ternary  exchange  systems,   the  following 
relation  must  hold  for  all  exchanger  phase  compositions 

dXi  =  0- 

This  criterion  holds  if 

P^^/^-  +  2i  P;,i/*  +  Zj  P^^/i    =  0,   for  n  =  1,   2,   3...  [54] 

where  P^^/J  is  the  n-th  parameter  in  the  regression  equation 
[52]   or   [53],   and  p^J/^c  and  P„*/i  are  the  corresponding  n-th 
parameters  in  the  ternary-based  regression  equations  having 
^^k(v)j/k    and  lnA:(v);c/i)  ,   respectively,  as  their  dependant 
variables . 

None  of  the  regression  parameter  estimates  deviated 
greatly  from  the  relation  expressed  in  equation   [54]  even 
though  the  restriction  was  not  included  in  the  calculation  of 
the  parameter  estimates.  Both  models  were  used  to  estimate 
the  exchanger  phase  composition.   For  both  the  cation  resin 
and  montmorillonite,   the  coefficient  of  determination  was 
greater  if  the  simpler  regression  model   (Equation   [52])  was 
used.  As  with  the  regression  equations  determined  from  binary 
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exchange  data,  the  parameter  estimates  were  adjusted  to 
adhere  to  [54].  The  adjusted  parameter  estimates  for  equation 
[52]   for  the  cation  exchanger  are  presented  in  Table  19  and 
for  montmorillonite  in  Table  20.  Since  the  parameter 
estimates  were  not  calculated  directly  from  least  squares 
regression  but  adjusted  to  adhere  to  equation  [54],  the 
coefficients  of  determination  from  the  original  polynomial 
regressions  are  not  presented. 

Calculation  of  Exchanger  Pha.gp  Activity  CoeffiriPnt 

The  exchanger  phase  activity  coefficients  were  estimated 
by  evaluating  the  expressions 

1 

In       =  — 7  (  (x/,  -  1)  lnJt„^,     .  +     Clnk  ^        dx^    )  r^si 

2j  ^  (V)i/j7  I  (V)i/j   Cl-^i  '    )r  [Zb] 


and. 


m  fj  =  jj  (x/,  ln;c,J,,/.  -  Jln;c,J,,/.  dx/,  )  [27] 

0 

where  In  k^^^^./j  was  estimated  by  the  quadratic  regression 

equation  calculated  from  the  binary  exchange  data.  Using  the 
calculated  exchanger  phase  activity  coefficients,  the 
exchanger  phase  activities  could  be  estimated  directly.  For 
binary  systems,   these  estimates  could  be  compared  to  the 
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Table  19.  Polynomial  model  parameter  estimates  from  least 
squares  regressions  of  the  natural  logarithm  of  adjusted 
Vanselow  selectivity  coefficients  on  cation  exchanger  phase 
equivalent  fractions  in  ternary  systems  with  AG  MP-50 
macroporous  cation  resin. 


ParamptPr 

Exchange 
selectivity 

Pot 

Pi 

P2 

1"  •'f(v)Mg/Ca 

-2 

.806 

3.256 

-0.236 

-3. 

213 

In  -?C(v)Na/Mg 

-0 

725 

2.135 

0.533 

-0. 

815 

In  >^(v)Ca/Na 

2  . 

128 

-3.763 

-0.415 

2. 

421 

Parameter  estimates  for  the  model 
^(v)i/j  =  Po  +  Pixi    +  p2^2  +  ^3x1X2,  where 


xi  -  exchanger  phase  Ca  equivalent  fraction 
X2  =  exchanger  phase  Na  equivalent  fraction. 


Table  20.  Polynomial  model  parameter  estimates  from  least 
squares  regressions  of  the  natural  logarithm  of  adjusted 
Vanselow  selectivity  coefficients  on  cation  exchanger  phase 
equivalent  fractions  in  ternary  systems  with  montm?rU?onUe 


Exchange 
selectivity 

Par;=)m(=>l-PT 

estimatpc; 

Pot 

Pl 

P2 

Pa 

^(v) Ba/NH4 

-2 

.363 

2.747 

-0. 024 

-0 

.551 

-^(v)  La/Ba 

2 

.372 

0.477 

0  .207 

-1 

2  92 

^ (v)NH4/La 

2 

.359 

-4.360 

-0.067 

1 

473 

re 


^(v)i/j  =  Po  +  Pixi    +  P2X2  +  P3X1X2,  whe 

xi  =  exchanger  phase  NH4  equivalent  fraction 
X2  -  exchanger  phase  Ba  equivalent  fraction 
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values  assumed  by  the  three  exchanger  phase  reference 
functions  considered  in  this  study;  the  mole  fraction,  the 
equivalent  fraction,   and  the  statistical  thermodynamic. 
Typical  plots  showing  the  relation  between  these  reference 
functions  and  the  calculated  exchanger  phase  activity  for  Na 
on  the  cation  resin  and  NH4  on  montmorillonite  are  presented 
in  Figures  10  and  11,   respectively.  The  closeness  of  the 
three  reference  functions  to  the  calculated  exchanger  phase 
activity  was  examined  by  calculating  the  square  root  of  the 
mean  square  error   (RMSE)   between  the  reference  function  and 
the  calculated  exchanger  phase  activity.  The  results  of  these 
calculations  are  presented  in  Table  21  for  the  cation  resin 
and  Table  22  for  montmorillonite. 

The  RMSE  statistic  was  used  extensively  in  this  study  to 
compare  the  relative  accuracy  of  predictions  using  chemical 
submodels  to  data  from  batch  and  column  experiments.  The  use 
of  statistics  to  validate  and  evaluate  computer  simulation 
models  does  not  appear  to  be  developed.  Given  the  paucity  of 
statistical  techniques,   its  was  decided  to  turn  to  regression 
analysis  for  guidance.   In  regression  analysis,  one  method  of 
validating  a  candidate  regression  model  is  to  collect  fresh 
data   (that  is  data  which  has  not  been  used  in  the  calculation 
of  the  parameter  estimates)   and  compare  the  r2  for  prediction 
^^prediction^   "^^^^  the  Coefficient  of  determination   (R2)for  the 
original  regression   (Montgomery  and  Peck,   1982) .   The  r2 

.        ,    ^  .        ,   ,  Prediction 

IS  defined  by 
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Table  25.  Square  root  mean  square  error  between  the  measured 
exchanger  phase  mole  fractions  on  AG  MP-50  macroporous 
cation  resin  and  exchanger  phase  mole  fractions  calculated  by 
using  polynomial  multiple  regression  fits  of  the  natural 
logarithm  of  adjusted  Vanselow  selectivity  coefficients  to 
exchanger  phase  equivalent  fractions. 


Data  used  in  rpgressinn 


Cation  predicted  Binary  Ternary 


— Root  mean 

square  error — 

Na 

0.080 

0.066 

Ca 

0.079 

0.065 

Mg 

0.063 

0.054 
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Table  26.   Square  root  mean  square  error  between  the  measured 
exchanger  phase  mole  fractions  on  montmorillonite  and 
exchanger  phase  mole  fractions  calculated  by  using  polynomial 
multiple  regression  fits  of  the  natural  logarithm  of  adjusted 
Vanselow  selectivity  coefficients  to  exchanger  phase 
equivalent  fractions. 


Data  used  in  reares.q  i  nq 


Cation  predicted 


Binary 


Ternary 


— Root  mean  square  error — 


NH4 


0.070 


0.042 


Ba 


0.011 


0.027 


La 


0.076 


0.065 
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11 


r2   .  ,     =       1  - 


Predictior  ~  [55. 

tfvi  -  y  f 


i=l 


where 

Yi  =  the  i-th  observation  of  the  dependant  variable 

y-hati  =  the  i-th  predicted  value  of  the  dependant 
variable 

y-bar    =  the  mean  of  all  observed  values  of  the 
dependant  variable. 

Since  the  total  sum  of  squares,   the  denominator  in 
equation   [56],  would  be  constant  for  all  models  evaluated 
with  fresh  data,   the  Rp'^^,,^^.^^  measure  would  be  inversely 
proportional  to  the  mean  square  error   (MSE  or  MSe) 

MSE     =  ^  

n 

For  qualitative  comparisons,   the  square  root  of  the  mean 
square  error   (RSME)   has  an  advantage  over  the  MSE,  because 
the  RMSE  measure  has  the  same  units  as  the  dependant 
variable.   Unlike  the  R^^,,^^^^^  measure,   both  the  MSE  and  RMSE 
do  not  include  information  on  the  total  variance  of  the 
dependant  variable.   It  should  be  recalled  that  the  total  sum 
of  squares  from  column  studies,   of  the  solute  concentrations 
in  the  column  effluent  concentrations  for  example,   would  be 
much  higher  than  the  prediction  error  sum  of  squares  (the 
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numerator  in  equation  [55]),   and  calculated  r2  could 

Prediction 

imply  a  unrealistic  level  of  accuracy  to  the  model  under 
evaluation.  All  of  the  statistics,  Rp^^^,^^.^^,  MSE  and  RMSE  are 
qualitative  measures  of  predictive  performance  and 
statistically  valid  test  of  hypotheses  cannot  be  calculated 
from  them.  Some  may  find  them,  however,   superior  to  plots  of 
measured  and  simulated  data  in  conveying  the  accuracy  of 
computer  simulations. 

No  exchanger  phase  reference  function  was  consistently 
closest  to  the  exchanger  phase  activity  in  binary  exchange  on 
the  cation  exchanger,   though  the  equivalent  fraction  failed 
to  give  the  best  estimate,   as  measured  by  RMSE,   of  any  of  the 
binary  systems.  For  all  of  the  binary  exchange  equilibria 
reported  by  Elprince  et  al .    (1980),  the  mole  fraction  was 
consistently  closer  to  the  calculated  exchanger  phase 
activity  and  predictions  due  to  the  equivalent  fraction 
reference  function  were  consistently  furthest.  Therefore,  for 
the  two  exchangers  considered  in  this  study,   the  mole 
fraction  reference  function  consistently  gave  a  close 
approximation  of  the  exchanger  phase  activity.     This  result 
suggested  that  the  mole  fraction  may  be  the  preferred 
reference  function  in  chemical  submodels  for  cation  transport 
models. 

The  ability  of  a  reference  function  to  approximate 
closely  the  exchanger  phase  activity  for  a  wide  range  of 
compositions  may  be  an  important  property  for  adherence  of 
the  corresponding  defining  relation  to  the  thermodynamic 
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theory  of  ion  exchange.  Practically,   and  more  importantly  for 
this  study,  the  preferred  reference  function  should  aid  in 
the  accurate  calculation  of  the  exchange  equilibria.  This 
ability  was  tested  for  the  three  reference  functions  and  for 
least  squares  regression  fits  of  the  In  k^^^^.^.  to  exchanger 
phase  equivalent  fractions  for  binary  and  ternary  exchange 
systems . 


Evaluation  of  Exchanger  PhasP  ReferPnr,^  Func1-innc, 

Program  BLCKBX  was  used  to  calculate  the  exchanger 
phase  mole  fractions  for  the  cation  resin  and  montmorillonite 
using  the  experimental  supernatant  solution  compositions,  the 
thermodynamic  exchange  constants  and  each  of  the  three 
exchanger  phase  reference  functions.  Comparisons  of  predicted 
and  measured  exchanger  phase  compositions  in  ternary  exchange 
are  presented  in  Figures  12  through  17  for  the  cation  resin 
and  montmorillonite. 

Several  of  the  plots  show  that  exchange  information  was 
not  gathered  in  the  for  a  full  range  of  exchanger  phase  mole 
fractions.  This  is  most  notable  for  Mg  on  the  cation  resin 
(Figure  14)   and  for  NH4  on  montmorillonite    (Figure  17).  with 
a  limited  number  of  equilibrations,   it  is  difficult  to  select 
correctly  the  equilibrating  solution  compositions  which  will 
give  exchange  information  for  the  range  of  exchanger  phase 
compositions,   especially  for  less  competitive  cations.  Both 
NH4  and  Mg  were  the  least  preferred  cations  in  their 
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respective  ternary  exchange  experiments  and  it  is  not 
surprising  that  exchange  data  was  lacking  for  higher 
exchanger  phase  mole  fractions  of  these  cations. 

As  before,  the  adequacies  of  various  reference 
functions  were  evaluated  by  calculating  the  RMSE  between 
predicted  and  measured  exchanger  phase  mole  fractions.  Tables 
23  and  25  show  these  RMSE  values  for  the  cation  resin  and 
montmorillonite,   respectively.  For  ternary  systems,  the  use 
of  the  equivalent  fraction  reference  function  appeared  to 
give  the  closest  estimate  of  Ca  and  Na  mole  fraction,  while 
none  of  the  RMSE  values  for  predicted  exchanger  phase  Mg  mole 
fractions  differed  greatly. 

In  the  binary  exchange  systems,  the  defining  relations 
for  the  homovalent  Ca-Mg  system  were  identical  and  there  were 
no  difference  between  the  predictions  due  to  the  three 
reference  functions.  The  equivalent  fraction  reference 
function  appeared  to  give  the  most  accurate  prediction  of 
exchanger  phase  composition  for  Ca-Na  systems,  while  use  of 
the  mole  fraction  reference  function  gave  the  best  estimate 
in  Mg-Na  exchange  systems. 

The  mole  fraction  reference  function,  which  gave  the 
best  approximation  of  the  calculated  exchanger  phase  activity 
for  exchange  on  montmorillonite  reported  by  Elprince  et  al. 
(1980),     also  gave  the  most  reliable  prediction  of  exchanger 
phase  composition  in  binary  and  ternary  exchange  experiments. 

The  exchanger  reference  function  which  yielded  the  most 
accurate  prediction  of  exchanger  phase  composition  was 
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different  for  the  two  experiments.  For  the  data  reported  by. 
Elprince  et  al.    (1978),  the  use  of  the  mole  fraction 
exchanger  phase  reference  function  consistently  gave  more 
accurate  estimates  of  exchanger  phase  compositions  as 
indicted  by  lower  RMSE  between  measured  and  predicted 
exchanger  phase  cation  mole  fractions.  For  exchange  with 
synthetic  exchange  resins,  the  exchanger  phase  composition 
predictions  using  equivalent  fraction  reference  function 
proved  to  be  the  most  reliable,  though  differences  among  the 
predictions  were  not  great.  While  it  would  be  tempting  to 
conclude  that  the  different  results  between  the  two 
experiments  were  due  to  the  nature  of  the  exchangers,  no 
conclusions  could  be  drawn  to  the  cause  of  the  different 
predictive  performance  of  the  reference  functions  on  the  two 
ion  exchangers  because  there  were  several  differences  between 
the  two  sets  of  experiments  aside  from  exchanger  type. 
Chloride  salts  were  used  by  Elprince  et  al.    (1978)  while 
perchlorate  salts  were  used  in  this  study.  The  system  studied 
by  Elprince  et  al.    (1978)   consisted  of  cations  with  valences 
of  1+,  2+,  and  3+,  while  the  cations  used  in  the  present 
study  had  valences  of  1+,  2+  and  2+.  The  equilibrating 
solutions  used  by  Elprince  et  al .  were  maintained  at  constant 
chloride  concentration  of  0.1  mol  dm-3,  while  the 
equilibrating  solutions  in  this  study  were  maintained  at  a 
constant  perchlorate  concentration  of  0.01  mol  dm-3.  Although 
it  may  be  suspected  that  the  two  exchanger  phase  reference 
function  performed  better  on  the  different  exchangers. 
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confounding  factors  eliminate  the  possibility  of 
demonstrating  conclusively  this  possibility.     The  prediction 
RMSE  using  the  equivalent  fraction  reference  function  were 
smallest  of  the  three  when  the  differences  between  the  three 
predictions  were  small.  The  differences  between  the 
predictions  were  greater  when  predicting  the  exchange 
equilibria  of  the  data  of  Elprince  et  al.    (1980).   In  that 
case,  the  predictions  made  using  the  mole  fraction  reference 
function  yielded  the  smallest  root  mean  square  error.  This 
trend  was  observed  in  evaluation  of  the  predictions  of 
cationic  solute  concentrations  in  the  column  effluent. 

Regression  equations  derived  from  binary  and  ternary 
data  were  evaluated  for  their  ability  to  predict  the  behavior 
of  the  ternary  systems.  The  tabulation  of  the  RMSE  between 
measured  exchanger  phase  mole  fractions  and  those  predicted 
by  the  use  of  variable  selectivity  coefficients  for  the 
cation  resin  and  montmorillonite  are  presented  in  Tables  24 
and  26,  respectively.  When  compared  to  the  best  of  suggested 
reference  functions,  the  use  of  the  binary-based  regression 
equations  reduced  the  root  mean  square  error  of  the 
prediction  by  about  half  when  the  resin  was  the  exchanger, 
but  much  less  with  the  montmorillonite  cation  exchanger.  With 
both  exchangers,  there  appeared  to  be  an  improvement  in  the 
root  mean  square  error  of  about  15%  by  using  ternary-based 
regression  equations  instead  of  the  binary-based  regression 
equations.  These  results  support,  though  weakly,  the 
contention  of  Chu  and  Sposito   (1980)  that  ternary  exchange 
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data  was  needed  to  describe  a  ternary  exchange  system.  It 
cannot  be  said  if  these  improvements  are  "significant" 
because  individuals,   researchers,  engineers,   or  planners 
would  have  to  judge  if  the  additional  effort  to  determine 
these  data  will  justify  the  expected  increase  in  accuracy. 

Regression  equations  which  relate  the  logarithm  of  the 
adjusted  Vanselow  selectivity  coefficient  to  the  exchanger 
phase  composition,  the  binary-based  regression  equations  and 
ternary-based  regression  equations,  have  simple  mathematical 
forms.  The  binary-based  regression  equations  were  degree  two 
polynomials  in  the  equivalent  fraction  of  one  of  the  two 
cations.  The  ternary-based  regression  equations  were  also 
degree  two  but  only  degree  one  in  the  equivalent  fraction  of 
either  of  the  two  cations  which  were  independent  variables. 
Attempts  to  include  higher  order  terms  in  the  regression 
analysis  did  not  improve  the  coefficient  of  determination  for 
these  regression  equations.  The  mathematical  simplicity  of 
the  ternary-based  regression  equations  suggests  that  their 
parameters  could  be  estimated  successfully  by  using 
binary-based  regression  equations. 

Binary-based  regression  equations  were  used  in  this 
study    calculate  the  exchanger  phase  activity  coefficients 
for  binary  exchange  using  equations   [26],    [27],   and  [51].  m 
principle,   regression  equations  derived  from  ternary  data 
could  be  used  to  calculate  exchanger  phase  activity 
coefficients  for  ternary  systems.  However,  the  derivation  of 
these  equations  are  beyond  the  scope  of  this  study.  Also 
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beyond  the  scope  of  this  study,  was  an  evaluation  of  the 
ability  of  the  regression  parameters  derived  from  binary  data 
to  estimate  the  parameters  of  the  ternary-based  regression 
equations.  If  exchanger  phase  activity  coefficients  could  be 
calculated  for  ternary  systems  as  a  function  of  the  adjusted 
Vanselow  selectivity  coefficient  and  if  the  ternary-based 
regression  coefficients  could  be  estimated  from  binary  data, 
then  equations  describing  the  exchanger  phase  activity 
coefficients  in  ternary  systems  could  be  derived  accurately 
from  binary  data. 


Column  Experiment? 


The  cation  concentrations  in  the  column  effluent 
solutions  are  presented  in  Tables  27  and  28.  As  noted  before. 
Run  A  consisted  of  a  saturation  of  the  column  with  0.005  mol 
dm-3    Ca(C104)2  to  which  61  cm3  of  a  0.0025  mol  dm-3 
Mg(C104)2,    0.005  mol  dm-3  NaC104  solution     was  pumped.  Run  B 
consisted  of  a  saturation  of  the  column  with  0.01  mol  dm-3 
NaC104  to  which  61  cm3  of  a  0.0025     mol  dm-3  Mg (0104)2,  0.0025 
mol  dm-3    ca (0104)2  solution  was  pumped.  The  parameters  used 
by  program  TABMODEL  to  simulate  these  column  experiments  are 
presented  in  Table  29. 
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Table  27.  Measured  column  effluent  mass  and  cation 
concentrations  in  Run  A;  which  consisted  of  a  pulse  of  a  0  05 
mol  dm  3  NaClO4-0.025  mol  dm-3  Mg(C104)2     solution  introducted 
into  a  column  saturated  with  a  0.05  mol  dm-3  Ca (0104)2 
solution . 


Effluent  Cation  concentration 
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mass  Na  Ca 
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Table  28.  Measured  column  effluent  mass  and  cation 
concentrations  in  Run  B;  which  consisted  of  a  pulse  of  a  0  25 
mol  dm-3  Ca(ClO4)2-0.025  mol  dm-3  Mg (0104)2  solution 
introducted  into  a  column  saturated  with  a  0.10  mol  dm-3 
NaC104  solution. 


Effluent  Cation  concentration 

fraction 
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Comparison  of  Tabular  Look-up  and  Valocrhi  Approache.q 

Two  criteria  were  used  to  evaluate  the  ability  of  the 
tabular  look-up  and  Valocchi  approaches  to  predict  the 
one-dimensional  transport  of  three  cations  through  porous 
media.  The  first  criterion  was  the  closeness  the  predicted 
first  moments  of  cation  concentration  were  to  those  of  the 
experimental  data.  The  first  moment  is  mean  time   (measured  in 
pore  volumes)   a  solute  pulse  remains  in  the  column.  The  first 
moment  is  defined  by 

Jvc  dv 

^'^  =  [55] 
Jc  dV 

where  c    is  the  solution  concentration  in  the  column 
effluent  and  V  is  the  pore  volume  value   (Villermaux, 1981) . 
The  second  criterion  was  the  RSME  between  measured  and 
predicted  cation  solution  concentrations  in  the  column 
effluent . 

The  measured  column  effluent  data  for  Run  A  are 
presented  in  Figure  18  along  with  the  predictions  of  an 
adaption  of  the  Valocchi  approach  and  predictions  using  the 
tabular  look-up  method  with  the  equivalent  fraction 
exchanger  phase  reference  function.  As  shown  in  Table  18,  the 
thermodynamic  exchange  constants  indicate  that  Ca  is 
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preferred  by  the  resin  over  Na  and  that  Ca  is  strongly 
preferred  by  the  resin  over  Mg.   The  value  calculated  for 
^^e^Mg/Na  was  0.524,   indicting  a  preference  by  the  cation  resin 
for  Na  over  Mg .   The  measured  Na  peak  was  approximately 
symmetric,   while  the  Mg  concentration  in  the  column  effluent 
had  the  shape  of  a  leading  peak,   an  indication  of  a  concave 
sorption  isotherm  faced  by  the  solute  in  contact  with  the 
porous  medium  (Schweich  and  Sardin,   1981) .  These  general  peak 
shapes  were  duplicated  successfully  by  each  of  the  modeling 
approaches . 

The  solute  pulse  first  moments  for  the  measured  data 
and  the  two  simulations  are  presented  in  Table  30.  Presented 
in  Table  31  are  root  mean  square  error  between  the  measured 
data  and  the  simulations  using  the  Valocchi  and  tabular 
look-up  methods.  As  can  be  seen  in  Figure  18,   there  were 
small  differences  between  the  two  predictions,   though  the 
results  of  the  tabular  look-up  approach  gave  more  accurate 
predictions  of  the  distribution  of  cation  concentrations  in 
the  effluent  as  measured  by  RMSE  than  did  the  Valocchi 
approach . 

The  calculation  of  exchange  equilibria  in  the  Valocchi 
approach  of  describing  cation  transport  through  porous  media 
differs  operationally  from  the  tabular  look-up  approach.  In 
the  Valocchi  approach,   the  concentrations  of  cations  in 
solution  determine  the  exchanger  phase  cation  molality,  while 
the  tabular  look-up  approach  uses  the  total  amount  of 
cations  in  both  solution  and  exchanger  phases  to  estimate  the 
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Table  31.   Square  root  mean  square  error  between  the  measured 
column  effluent  solution  cation  equivalent  fraction  and  those 
predicted  using  Valocchi  and  tabular  look-up  approaches  in 
the  development  of  cation  transport  models 


Simulation  approach  uf^f-<^ 


ValQCChi   Tabular  look-up 

Exchanger  phase 

^  ^  .       ,  reference  funrtiop  n^^.^ 

cation  in 

.  Mole  Equivalent 


effluent  fraction 


fraction 


-RMSE+   (Solution  equivalent  fraction) - 

Ca  n-?^?  ^-^'^  0.073 

t  0-138  0.187 

^5  0.143  0.132  0.151 


^Square  root  mean  square  error. 
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equilibrium  distribution  of  the  cations  between  the  two 
phases.  The  latter  approach  has  greater  conceptual  appeal, 
but  judging  from  the  results  of  these  column  simulations,  the 
predictive  improvement  in  abandoning  the  Valocchi  method  for 
the  tabular  look-up  approach,  while  real,  was  modest. 
Further  experiments  with  exchangers  with  different 
selectivities  and  different  experimental  initial  and  boundary 
conditions  from  the  present  study  would  be  needed  to 
determine  the  practical  implications  of  these  differences  in 
operational  computation  of  cation  exchange. 

Evglu^tion  Of  Exchanger  Pha.<^e  Reference  Functinn.^ 

One  Objective  of  this  study  was  to  determine  which  of 
several  exchanger  phase  reference  functions  provided  the  most 
accurate  description  of  multicomponent  cationic  solute 
movement  in  porous  media  when  used  in  a  numerical  solute 
transport  model.  As  was  discussed  earlier,  reference 
functions  and  selectivities  chosen  for  use  in  solute 
transport  models  assume,  usually  implicitly,  that  the 
corresponding  reference  functions  approximate  solute 
exchanger  phase  activity.  From  a  practical  viewpoint, 
however,   the  more  important  property  of  a  particular 
exchanger  phase  reference  function  is  the  predictive  accuracy 
obtained  when  it  is  incorporated  in  a  solute  transport  model. 
If  a  specific  exchanger  phase  reference  function  gives  a  more 
accurate  estimate  of  the  exchanger  phase  activity  than  other 
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reference  functions,   such  an  exchanger  phase  reference 
function  would  be  expected  to  provide  more  reliable 
simulation  results  when  used  in  a  solute  transport  model. 

Two  sets  of  column  data  were  used  in  the  evaluation  of 
the  three  exchanger  phase  reference  functions  used  in  this 
study.  The  first  column  data  were  those  generated  in  Run  A  of 
this  study.  The  second  were  those  presented  by  Lai  et  al. 
(1978)   for  a  column  packed  with  Yolo  loam  soil.  While  they 
used  a  column  packed  with  soil  rather  than  a  quartz  sand- 
cation  resin  mixture,  the  experimental  boundary  and  initial 
conditions  investigated  by  Lai  et  al.    (1978)  were  similar  to 
Run  A.  Lai  et  al.    (1978)   equilibrated  their  soil  column 
packing  to  0.125  mol  Ca (CH3COO) 2  dm-3    solution  and  introduced 
a  0,55  pore  volume  pulse  of  a  0.0625  mol  MgCl2-0,125  mol 
NaCl  dm"-^  solution. 

The  measured  Run  A  effluent  data  and  those  predicted  by 
the  use  of  the  mole  fraction  reference  function  are  presented 
in  Figure  19,  These  data  are  presented  again  with  the 
prediction  using  the  statistical  thermodynamic  the  equivalent 
fraction  reference  function  in  Figure  20  and  with  the 
prediction  using  the  equivalent  fraction     reference  function 
in  Figure  21.  The  differences  between  the  effluent 
predictions  and  the  predictions,   as  measured  by  the  root  mean 
squared  error,   are  presented  in  Table  32.  While  differences 
among  the  predictions  were  not  great,   the  statistical 
thermodynamic  reference  function  gave  the  most  accurate 
predictions  for  Na  and  Ca  concentrations  in  the 
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Table  34.  Calculated  effluent  cation  solution  mass-balance 
errors  for  five  implementations  of  tabular  look-up  approach 
to  the     simulation  of  Case  A. 


Implementation  type  Total  Proportion  of 

mass-balance     total  error  for 
error  each  cation 


Cation 


Na  Ca  Mg 


Mole  fraction 
Equivalent  fraction 
Statistical  thermodynamic 
Binary  regression 
Ternary  regression 


-percent-    — percent  of  total — 


7 

.36 

7.25 

3. 

20 

89 

.55 

1 

.13 

19.26 

1 . 

82 

78 

.  93 

1 

.02 

3.37 

1. 

43 

95 

.20 

7 

.73 

6.58 

3. 

35 

90 

.07 

0 

.78 

1 .79 

9. 

06 

89 

.15 
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columneffluent  while  the  equivalent  fraction  reference 
function  gave  the  most  accurate  predictions  for  Mg 
concentrations  in  the  column  effluent.  The  measured  and 
predicted  concentration  first  moments  are  presented  in  Table 
33.  None  of  the  reference  functions  gave  a  superior 
prediction  of  the  first  moment  of  the  cation  concentration  in 
the  column  effluent.  The  mass-balance  errors  in  the  computer 
simulations  based  on  the  various  tested  models  are  presented 
in  Table  34. 

The  measured  effluent  data  presented  by  Lai  et  al. 
(1978)  and  those  predicted  by  the  use  of  the  mole  fraction 
reference  function  are  presented  in  Figure  22.  The  Na  peak 
has  a  slight  leading  shape,  indicating  that  Na  was  less 
preferred  by  the  soil  than  the  other  two  cations,  while  the 
Mg  peak  had  a  trailing  shape  indicating  that  Mg  was  preferred 
by  the  soil  over  its  counterions  in  the  system.  The  column 
effluent  data  are  presented  again  with  the  prediction  using 
the  equivalent  fraction  reference  function  in  Figure  23  and 
with  the  prediction  using  the  statistical  thermodynamic 
reference  function  in  Figure  24.  While  the  predictions  using 
the  mole  fraction  and  equivalent  fraction  reference  functions 
were  roughly  equivalent  and  closely  match  the  measured  data, 
the  equivalent  fraction  predictions  were  different  from  the 
others  and  predicted  the  column  response  much  more  poorly 
that  those  of  the  other  two  reference  functions.  This 
difference  in  predictive  ability  was  most  noticeable  for  the 
Na  peak.  Unlike  the  data  developed  from  the  cation  resin  used 
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in  this  study,  the  Na  and  Mg  peaks  presented  by  Lai  et  al. 

(1978)   for  soil  appeared  to  have  separated 
chromatographically  soon  after  being  introduced  into 
thecolumn.  Since  the  Na  peak  represents  essentially  a  binary 
heterovalent  system,   it  is  suggested  that  the  equivalent 
fraction  reference  function  may  be  less  accurate  for 
predicting  the  transport  of  cations  in  binary  heterovalent 
systems . 

In  simulations  of  the  column  effluent  data  reported  by 
Lai  et  al.    (1978),  the  use  of  the  mole  fraction  and 
statistical  thermodynamic  reference  functions  gave  excellent 
descriptions  of  their  data.  The  use  of  the  equivalent 
fraction  reference  function  resulted  in  substantial  errors  in 
effluent  cation  concentration  predictions.  As  with  the 
predictions  of  batch  cation  exchange  equilibria,   in  some 
experimental  conditions  all  three  of  the  reference  functions 
tested  here  predicted  column  effluent  cation  concentration 
distributions  with  equivalent  accuracy.  In  other  conditions, 
however,  the  mole  fraction  was  superior  to  the  equivalent 
fraction  reference  function.  It  is  suggested,  therefore,  that 
while  the  mole  fraction  reference  function  did  not  always 
give  the  most  accurate  predictions  of  exchange  equilibria,  it 
appeared  to  be  the  more  reliable  by  not  giving  large  errors 
under  certain  conditions.  The  use  of  the  mole  fraction 
exchanger  phase  reference  function  failed  to  be  the  most 
accurate  in  situations  where  modest  differences  existed 
between  predictions  made  by  using  the  various  exchanger  phase 
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reference  functions.   In  those  cases  where  substantial 
differences  existed,  the  mole  fraction  was  superior  in  the 
cases  evaluated  by  this  study.  An  exception  to  this 
observation  is  in  the  column  effluent  data  of  Lai  et  al. 
(1978),  where  predictions  using  the  statistical  thermodynamic 
reference  function  were  roughly  as  accurate  as  those  using 
the  mole  fraction  reference  function. 

Effect  Of  Variable  .Se1ectiv.1tv  CoeffiriPn^s  on  PrPHirt^n 

Column  Respon.qp 

All  solute  transport  models,   indeed  all  mathematical 
models,  have  some  error  in  their  predictions;  models  of 
solute  cations  in  heterovalent ,  multicomponent  systems 
require  complex  chemical  submodels  and  may  have  more  error 
than  solute  transport  models  describing  systems  with  more 
simpler  adsorption  behavior.  A  problem  addressed  in  this 
study  was  the  degree  to  which  improvements  in  the  chemical 
submodel  should  improve  overall  model  predictions  of  solute 
transport.  For  example,   assume  that  by  using  variable  rather 
than  fixed  exchange  selectivities  as  input  parameters  in  the 
model,  the  error  in  the  predictions  of  batch  exchange 
equilibria  was  decreased  by  50%.  One  would  not  necessarily 
know,  however,   if  this  improvement  in  accuracy  would  be 
duplicated  if  variable  exchange  selectivities  replaced  fixed 
exchange  values  in  a  cation  transport  model. 
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To  this  end,  regression  equations  which  related  the 
natural  logarithm  of  the  adjusted  Vanselow  selectivity 
coefficient  to  exchanger  phase  composition  were  used  to 
predict  the  exchanger  phase  composition  determined  in  the 
batch  experiments.  Two  sets  of  regression  equations  were 
calculated;  one  set  was  based  on  binary  exchange  data  while 
the  second  set  was  based  on  the  results  of  ternary  exchange 
experiments.  These  sets  of  equations  were  referred  to  as  the 
"binary-based  regression  equations"  and  the  "ternary-based 
regression  equations",   respectively.  Chu  and  Sposito  (1980) 
argued  forcefully  that  ternary  ion  exchange  data  were 
necessary  to  accurately  describe  ternary  exchange  equilibria. 
Increased  accuracy  of  these  two  types  of  regression  equations 
in  both  batch  and  solute  transport  systems  should  indicate 
the  predictive  advantage  from  the  use  of  ternary  exchange 
data. 

The  measured  effluent  data  from  Run  A,  those  predicted 
by  the  use  of  the  mole  fraction  reference  function  with 
constant  selectivity  coefficients,  and  those  predicted  by 
using  binary-based  regression  equations  are  presented  in 
Figure  25.  The  measured  effluent  data,  those  predicted  by  the 
use  of  the  mole  fraction  reference  function  with  constant 
selectivity  coefficients,  and  those  predicted  by  selectivity 
coefficients  predicted  by  regression  equations  obtained  by 
fits  to  ternary  data  are  presented  in  Figure  26.  The 
predictions  due  to  binary  and  ternary  regressions  are 
presented  in  Figure  27.  The  RMSE  of  these  predictions  are 
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presented  in  Table  32,   and  the  cation  column  effluent 
concentration  distribution  first  moments  are  presented  in 
Table  33.  Predictions  of  column  effluent  cationic  solute 
concentrations  were  improved  by  as  much  as  half  by  using 
binary-based  regression  equations.  No  improvement  was  found 
by  using  ternary-based  regression  equations.  The  predictions 
using  these  regression  equations  were,   if  anything,  worse 
than  those  using  binary-based  regression  equations.  This 
result  suggested  that  while  there  was  an  improvement  in  the 
prediction  of  ternary  batch  exchange  equilibria  by  using 
ternary  exchange  data,  there  was  no  similar  improvement  in 
the  prediction  of  solute  transport  by  incorporating  the 
ternary-based  regression  equations  in  chemical  submodels  of 
solute  transport  models. 

SuaaestPd  Future  Rp.qp;qr-r->^ 

Results  of  this  study  suggest  several  avenues  for 
further  research.  First,  under  some  experimental  conditions 
the  use  of  the  mole  fraction  exchanger  phase  reference 
function  gave  predictions  superior  to  the  other  tested 
reference  functions.  This  study,  however,   did  not  clarify 
under  what  specific  conditions  the  predictions  following  from 
the  various  reference  functions  will  give  similar  predictions 
of  exchange  equilibria  and  under  what  conditions  they  will 
give  significantly  different  predictions  in  batch  and  column 
experiments . 
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Second,  ternary  exchange  information  in  this  work 
conveyed  little  predictive  improvement  to  ternary  batch  and 
no  improvement  to  predictions  of  a  column  experiment.  Further 
research  should  address  the  possibility  that  binary  exchange 
data  is  adequate  for  the  description  of  quaternary  and  higher 
order  multication  systems. 

Third,  the  use  of  thermodynamic  description  of  ion 
exchange  reactions  as  chemical  submodels  in  solute  transport 
models  assumed  implicitly  that  the  kinetics  of  the  exchange 
process  did  not  affect  significantly  the  transport  of  the 
cations.  Kinetic  effects  should  be  most  noticeable  at  high 
interstitial  velocities  as  were  used  in  these  column 
experiments.  This  study  has  shown  that  accurate  descriptions 
of  column  response  could  be  made  using  only  the  thermodynamic 
description  of  cation  exchange,  even  at  flow  rates  greater 
than  typically  imposed  on  soil  column  studies.  Still,  this 
study  has  not  shown  if  the  model  predictions  could  be 
improved  by  including  kinetic  effects  in  the  chemical 
submodel.  Resolution  of  this  possibility  is  needed  in  future 
experiments . 


APPENDIX  A 
LIST  OF  SYMBOLS  USED  IN  THE  TEXT 


Note:  The  current  recommendations  of  the  International  Union 
of  Pure  and  Applied  Chemistry  (lUPAC)  was  used  as  the 
guide  for  the  symbols  and  nomenclature  used  in  this 
dissertation   (MacGlashen,   1979;  Irving,   1972).  The  only 
manor  deviations  from  the  lUPAC  recommendations  is  the 
use  of  the  apostrophe  after  a  subscript  to  indicate  that 
the  unit  of  measurement  is  on  an  equivalent  basis  (e.q. 


e 
a,- 


s 

c 


s 


s 


D 


Cj^  and  c_.,  are  the  solution  concentrations  of  "i"  in 

units  of  mol  dm-3    and        mol  dm-3,   respectively)  ,  the 

use  of  the  superscript  """e"  to  refer  to  the  exchanger 
(solid)  phase,  and  the  use  of  the  superscript  "s"  to 
refer  to  the  solution   (liquid)  phase. 

=  activity  of  solute  i  in  the  exchanger  phase 
s 

=  activity  of  solute  i  in  the  solution  phase 

C        =  cation  exchange  capacity  of  the  porous  medium 
(mol(  +  )  kg-i) 


i        -  ion  i  solution  concentration   (mol  dm-3) 

=  ion  i  solution  equivalent  concentration   (~  mol  dm-3) 

"7  '  ' 


Cv,  -  sum  of  the  charge  born  by  all  cations  in  solution 
(mol(  +  )  dm-3) 


-  Coefficient  of  longitudinal  dispersion   (dm2  s"!) 

=  Denominator  appropriate  for  an  exchanger  phase 
reference  function 
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di        =  distribution  coefficient   (dm^  kg"l) 
e 

d,        =  - 

e  =  Superscript  designating  the  exchanger  phase 

feq(i)=  equivalence  factor  of  cation  i 

fi        =  exchanger  phase  activity  coefficient 

g         =  Gas  phase 

1         =  a  designation  for  cation  i 

Ic  =  ionic  strength  of  the  electrolyte  solution  (concentration 
basis) 

j  =  a  designation  for  cation  j 

^(Rj)i/j  "     Rubin- James  selectivity  coefficient 

^(V)i/j  "    Vanselow  selectivity  coefficient 

^(VAL)  i/j  "    Valocchi  selectivity  coefficient 
,  a  _ 

^(GT)i/j  -  adjusted  Gaines-Thomas  selectivity  coefficient 
,  a  _ 

^{ST)i/j  -  adjusted  Statistical  Thermodynamic  selectivity 
coefficient 

^(V)i/j  -  adjusted  Vanselow  selectivity  coefficient 
,  a  _ 

i/j     ~        adjusted  exchange  selectiviy  coefficient 

^Cj      =  concentration  equilibrium  constant  for  the  production 
of  derived  species  j 

Kd        =  linear  partition  coefficient 

^exi/j  =  thermodynamic  exchange  constant 


=  an  exchange  selectivity  coefficient 
=  differential  operator 
52  5 

=  zero-order  rate  coefficient  of  production  in  solid (s 

=  zero-order  rate  coefficient  of  production  in 
solution   (mol  dm"3  s-l) 

=  first-order  rate  coefficient  of  decay  in  solid  (s"!) 

=  first-order  rate  coefficient  of  decay  in  solution 
(s-l) 

=  exchanger  phase  concentration,  mol  kg-l 

=  ion  i  solid  solution  equivalent  molality   {—  mol  kg-i 
=  number  of  cations  in  the  exchange  system 
=  number  of  moles  of  i  in  the  solid  phase 

=  number  of  moles  of  i  in  the  liquid  phase 


e 


e  ,  e 
ZiDj^  +  zj  rij 


p  b     +  e  cl 


-  total  concentration  of  cation  i  in  exchanger  and 
solution  phases   (mol  dm-3) 

=  number  of  moles  of  solute  i  in  the  solution  phase 
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R         =  molar  gas  constant   (J  mol"i) 
r(^)     =  a  reference  function 

e 

=  exchanger  phase  reference  function 

s 

rj_        =  solution  phase  reference  function 

s         =  Solution  phase 
T         =  temperature  (K) 
t         =  time  (s) 

u         =  interstitial  velocity   (dm  s'^) 

w  =  subscript  for  the  solute 

e  _ 

=  mole  fraction  of  cation  i  in  the  exchanger  phase 

e 

x^,  =  equivalent  fraction  of  cation  i  in  the  exchanger  phase 

Yi  =  multiplicative  factor  suggested  by  Krishnamoorthy  et 
al.    (1948)   equal  to   (Zi  +  l)/2 

z  =  length  in  the  z  direction  (dm) 

Zi  =  the  absolute  value  of  the  valence  of  ion  i 

Zj  =  the  absolute  value  of  the  valence  of  ion  j 

Zk  =  the  absolute  value  of  the  valence  of  ion  k 

At  =  finite  difference  of  time 

A;<  =  finite  difference  of  variable  x 

Ay  =  finite  difference  of  variable  y 

Az  =  finite  difference  of  variable  z 

Yi  =  solution  phase  activity  coefficient 

M'  =  chemical  potential 
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e 

-  exchanger  phase  cation  i  chemical  potential 

eo 

\i  I       -  standard  state  exchanger  phase  cation  i  chemical 
potential 
eo  _ 

Ji„       -  standard  state  exchanger  phase  solvent  chemical 
potential 

e 

=  exchanger  phase  solvent  chemical  potential 
s  _ 

-  solution  phase  cation  i  chemical  potential 

So 

M  j;       -  Standard  state  solution  phase  cation  i  chemical 
potential 

0         =  volumetric  water  content   (dm^  dm"^) 
Pb       =  porous  medium  bulk  density   (kg  dm-3) 
C*        =  the  reference  state 


APPENDIX  B 


LISTING  OF  COMPUTER  PROGRAM 


PROGRAM  TBLMKR 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK 
SCCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI  IRF 
&ISBMKR,    ITRLIM,    JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,   RKKI,   RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 • 101 ) 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

CHARACTER  FINPUT*50, FOUTPT*50 

WRITE (*,*)    'INPUT  FILE  NAME' 

READ {*, 140)  FINPUT 

OPEN (UNIT=10,  FILE=FINPUT, STATUS='OLD ' ) 
WRITE (*,*)    'OUTPUT  FILE  NAME' 
READ (*, 140)  FOUTPT 

OPEN (UNIT=4 ,  FILE=FOUTPT, STATUS= ' NEW ' ) 

C 

c 

0 
0 

c 
c 

C  THE  SUBROUTINES 

C 

c 

c 
c 

c 

C  BLCKBX  =  THIS  SUBROUTINE  CALCULATES  THE  PARTITIONING  OF  THREE  CATIONS 
C  BETWEEN  THE  EXCHANGER  AND  SOLUTION  PHASES.  ^^^lU^Jb 

C  CCALC  =  USING  THE  PROPORTION  OF  CEC*BULKD  AND  THETA*TOTNOR  THIS 

C  SUBROUTINE  GIVES  A  FIRST  ESTIMATE  OF  THE  DISTRIBUTION  OF 

C  THREE  IONS  BETWEEN  THE  TWO  PHASES 

C^DENOM  =    CALCULATES  THE  DENOMINATOR  OF  THE  REFERENCE  FUNCTION  CHOSEN 

C  APPROXIMATE  THE  EXCHANGER  PHASE  ACTIVITY. 

C  F  =  CALCULATES  THE  EXCHANGER  PHASE  ACTIVITY  COEFFICIENTS 

C  FCALC  =  CALCULATES  THE  EXCHANGER  PHASE  MOLE  FRACTION  FOR  EACH  CATION 
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C  GAMMA 

C 

C 

C  INPUT 
C 

C  I THEN J 
C 
C 
C 
C 
C 
C 
C 


CALCULATES  THE  SOLUTION  PHASE  ACTIVITY  COEFFICIENTS  FOR 
EACH  CATION 

INPUTS  PARAMETER  FILE 


(ALSO  ITHENK,   JTHENI,   JTHENK,   KTHENI,   KTHENJ)  THESE 
SUBROUTINES  USE  THE  DEFINING  RELATION  OF  THE  VANSELOW 
SELECTIVITY  COEFFICIENT  TO  ESTIMATE  ONE  OF  THE  CATIONS 
THE  SECOND  IS  ESTIMATED  BY  DIFFERENCE   (THE  THIRD  IS  FIXED) 
IF  THE  ESTIMATE  OF  THE  SECOND  ION  IS  NEGATIVE,   THEN  THE 
ION  CONCENTRATION  VALUES  REVERT  TO  THEIR  PREVIOUS  VALUES 
AND  ANOTHER  SUBROUTINE  IS  TRIED.   IF  THE  ESTIMATE  OF  THE 
SECOND  ION  IS  POSITIVE.   THE  TWO  IONS  ARE  ESTIMATED 
ITERATIVELY  ^i^±r.u 

C  UNTIL  A  CONVERGENT  SOLUTION  IS  REACHED 

C 

C  KCALC    =  CALCULATES  THE  VANSELOW  SELECTIVITY  COEFFICIENT   (ADJUSTED  FOR 

C  SOLUTION  PHASE  ACTIVITY  COEFFICIENTS) 

C 

C  OUTPUT  =  OUTPUTS  RESULTS,   NOT  CALLED  IN  CURRENT  IMPLIMENTATION 
0 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


THE  VARIABLES 


C 

c 

C  BULKD 
C 

C  CBARI 
C 

C  CBARJ 
C 

C  CBARK 
C 

C  CCI 
C 

C  CCJ 

c 

C  CCK 
C 

C  CCTN 


THE  VARIABLES  IN  COMMON  BLOCK  STEVES 
THE  BULK  DENSITY  (KG/DECIMETER'^ 3) 
EXCHANGER  PHASE  CONCENTRATION  (MOL/KG) 
EXCHANGER  PHASE  CONCENTRATION  (MOL/KG) 
EXCHANGER  PHASE  CONCENTRATION  (MOL/KG) 

THE  TOTAL  AMOUNT  OF  ION  I  IN  THE  SYSTEM  (MOL/DECIMETER^S) 
THE  TOTAL  AMOUNT  OF  ION  J  IN  THE  SYSTEM  (MOL/DECIMETER-3) 
THE  TOTAL  AMOUNT  OF  ION  K  IN  THE  SYSTEM  (MOL/DECIMETER^S) 
THE  SUM  OF  ALL  CATIONIC  CHARGE  IN  THE  SYSTEM 


(MOL  (+)  /DECIMETER'^ 3) 


C  CECFIX  =  THE  CATION  EXCHANGE  CAPACITY (MOL (+) /KG) 

C  CI  =  SOLUTION  PHASE  CONCENTRATION  (M0L/DECIMETER"3) 
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C 

C  CJ  =  SOLUTION  PHASE  CONCENTRATION  (MOL/DECIMETER'^3) 

C  CK  =  SOLUTION  PHASE  CONCENTRATION  (M0L/DECIMETER-"3) 

C  CTOTAL  =  THE  TOTAL  SOLUTION  NORMALITY (MOL (+) /DECIMETER^S) 

C  FI  =  THE  EXCHANGER  PHASE  ACTIVITY  COEFFICIENT  FOR  ION  I 

(DIMENSIONLESS) 
C 

C  FJ  =  THE  EXCHANGER  PHASE  ACTIVITY  COEFFICIENT  FOR  ION  J 

(DIMENSIONLESS) 
C 

C  FK  =  THE  EXCHANGER  PHASE  ACTIVITY  COEFFICIENT  FOR  ION  K 

(DIMENSIONLESS) 
C 

C  EPFRCI  =  THE  ION  I  EXCHANGER  PHASE  MOLE  FRACTION  (DIMENSIONLESS) 

C  EPFRCJ  =  THE  ION  J  EXCHANGER  PHASE  MOLE  FRACTION  (DIMENSIONLESS) 

C  EPFRCK  =  THE  ION  K  EXCHANGER  PHASE  MOLE  FRACTION  (DIMENSIONLESS) 

C  EXPFRC  =  PROPORTION  OF  THE  SYSTEM  TOTAL  CATIONIC  CHARGE  FOUND  IN  THE 
C  EXCHANGER  PHASE  (DIMENSIONLESS) 

C  GAMMAI  =  THE  SOLUTION  PHASE  ACTIVITY  COEFFICIENT  FOR  ION  I 
(DECIMETERS 3 /MOL) 
C 

(D??^?ER'^3/^Lr"^'°''  COEFFICIENT  FOR  ION  J 

C 

(DES?E;>^3^Lr'^'°''  ^"^^^  COEFFICIENT  FOR  ION  K 

C 

C  IDEBUG  =  MARKER  CONTROLLING  DEBUGGING  PRINTOUT 
C  IDEBUG  LT  1     NO  DEBUGGING  PRINTOUT 

c 

C  IJ  =  FLAG  WHICH  INDICTES  IF  THE  BINARY  EQUILIBRIUM 

r  ^T^^^  ^^^^  I  AND  J  (K  FIXED)   HAS  BEEN  SOLVED  FOR 

C  THE  CURRENT  ITERATION 

C 

§  =  S^frSo^L^"'^  '^^'^  °^  ^^^^  EQUILIBRIA  SOLVED  FOR 
C  2  ™^  ^TE^TION  IS  COMPLETE  WHEN 

C 

C  =  ITl^^™^  "^"^^^  REFERENCE  FUNCTION  HAS  BEEN 

^  CHOSEN  TO  REPRESENT  THE  CATION  EXCHANGER  PHASE  ACTIVITY 

C  MOLE  FRACTION  =  1;   EQUIVALENT  FRACTION  =  2- 

C  STATISTICAL  THERMODYNAMICAL  =  3 

C  ISBMRK  =  MARKER  INDICATING  WHICH  OF  THE  ATHENB  IS 
C  CALLING  THE  DENOM  SUBROUTINE 

C  ITHENJ  =  1;   ITHENK  =  2;   JTHENI  =  3- 

C  JTHENK  =  4;  KTHENI  =  5;  KTHENJ  =  6 


IDEBUG  GE  1     DUBUGGING  STATEMENTS  WRITTEN 
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SOLUTION  PHASE  ACTIVITY  COEFFICIENTS, (DIMENSIONLESS) 


C  ITRLIM  =  THE  MAXIMUM  NUMBER  OF  ITERATIONS  ALLOWED  IN  ANY 

C  SUBROUTINE  OR  THE  MAIN  PROGRAM 

C 

C  JK  =  FLAG  WHICH  INDICTES  IF  THE  BINARY  EQUILIBRIUM 

C  BETWEEN  IONS  J  AND  K   (I  FIXED)   HAS  BEEN  SOLVED  FOR 

C  THE  CURRENT  ITERATION 

C 

C  KI  =  FLAG  WHICH  INDICTES  IF  THE  BINARY  EQUILIBRIUM 

C  BETWEEN  IONS  K  AND  I    (J  FIXED)   HAS  BEEN  SOLVED  FOR 

C  THE  CURRENT  ITERATION 
C 

C  RKIJ  =  THE  KVIJ  VANSELOW  SELECTIVITY  COEFFICIENT  ADJUSTED  TO 
INCORPORATE 

C  SOLUTION  PHASE  ACTIVITY  COEFFICIENTS, (DIMENSIONLESS) 

C  RKIJO  =  THE  THERMODYNAMIC  EXCHANGE  CONSTANT,   KIJ  (DIMENSIONLESS) 

C  RKJK      =  THE  KVJK  VANSELOW  SELECTIVITY  COEFFICIENT  ADJUSTED  TO 
INCORPORATE 
C 
C 

C  RKJKO     =  THE  THERMODYNAMIC  EXCHANGE  CONSTANT,   KJK  (DIMENSIONLESS) 

C  RKKI      =  THE  KVJK  VANSELOW  SELECTIVITY  COEFFICIENT  ADJUSTED  TO 
INCORPORATE 

C  SOLUTION  PHASE  ACTIVITY  COEFFICIENTS,  (DIMENSIONLESS) 

C  RKKIO     =  THE  THERMODYNAMIC  EXCHANGE  CONSTANT,   KKI  (DIMENSIONLESS) 

C  SNPFRC  =  PROPROTION  OF  THE  TOTAL  CATIONIC  CHARGE  FOUND  IN  THE  SOLUTION 

C  PHASE  (DIMENSIONLESS) 

C 

C  TABLE    =  ARRAY  OF  SOLUTION  PHASE  COMPOSITIONS  (DIMENSIONLESS) 

C  TFRACI  =  PROPROTION  OF  THE  TOTAL  SYSTEM  CATIONIC  CHARGE 
C  BORN  BY  THE  I  IONS  (DIMENSIONLESS) 

C  TFRACJ  =  PROPROTION  OF  THE  TOTAL  SYSTEM  CATIONIC  CHARGE 

C  BORN  BY  THE  J  IONS  (DIMENSIONLESS) 

C 

C  TFRACK  =  PROPROTION  OF  THE  TOTAL  SYSTEM  CATIONIC  CHARGE 
C  BORN  BY  THE  K  IONS  (DIMENSIONLESS) 

C  THETA    =  THE  VOLUMETRIC  WATER  CONTENT (DECIMETER- 3 /DECIMETERS 3) 

C  TOL        =  ITERATIVE  TOLERENCE  FOR  ROOT  MEAN  SQUARE  ERROR 
C  OF  THE  EXCHANGER  PHASE  ESTIMATES 

c 

C  DENOM  =    THE  DENOMINATOR  IN  THE  REFERENCE  FUNCTION  FOR  EXCHANGER  PHASE 

r  SSt7w3./"^^^^  VANSELOW:MOL/KG,  GAINES-THOMAS: 

MOL(  +  )/KG) 

C 

C  ZI  =  THE  VALENCE  OF  ION  I  (DIMENSIONLESS) 

C  ZJ  =  THE  VALENCE  OF  ION  J  (DIMENSIONLESS) 

C  ZK  =  THE  VALENCE  OF  ION  K  (DIMENSIONLESS) 
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C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


MAIN  PROGRAM 


N 

RN 

I 

ILOWER  = 
lUPPER  = 


THE  DIMENSION  OF  THE  TABLE  TO  BE  WRITTEN  TO  ARRAY  "TABLE" 
RATIONAL  VALUE  OF  THE  INTEGER  VARIABLE  "N" 
OUTER  DO  LOOP  VARIABLE 

DEBUGGING  AID,  SUBROUTINE  BLCKBX  WILL  NOT  BE  CALLED  IF  THE 
LOOP  VARIABLE  I  IS  LESS  THAN  THIS  VARIABLE 

DEBUGGING  AID.  SUBROUTINE  BLCKBX  WILL  NOT  BE  CALLED  IF  THE 
LOOP  VARIABLE  I  IS  LESS  THAN  THIS  VARIABLE 


J 

JLOWER 


=  INNER  DO  LOOP  VARIABLE 

=  DEBUGGING  AID.   SUBROUTINE  BLCKBX  WILL  NOT  BE  CALLED  IF  THE 
LOOP  VARIABLE  J  IS  LESS  THAN  THIS  VARIABLE 


JUPPER  = 


NJ 


DEBUGGING  AID.   SUBROUTINE  BLCKBX  WILL  NOT  BE  CALLED  IF  THE 
LOOP  VARIABLE  I  IS  GREATER  THAN  THIS  VARIABLE 

N  MINUS  J,   USED  IN  CONSTRUCTION  OF  UPPER  DIAGONAL  OF  "TABLE" 


SUBROUTINE  BLCKBX  AND  THE  ATHENB  SUBROUTINES 


ITERl 
ITER2 
CBARIO 
CBARJO 
CBARKO 
CBARIl 
CBARJl 
CBARKl 
CBARI2 
CBARJ2 
CBARK2 
CBICI 

CIO 


=  THE  COUNTER  FOR  THE  INNER  ITERATIVE  LOOP 
=  THE  COUNTER  FOR  THE  OUTER  ITERATIVE  LOOP 

=  HOLDS  ORIGINAL  VALUE  OF  CBARI  DURING  ITERATIVE  SOLUTION 

=  HOLDS  ORIGINAL  VALUE  OF  CBARJ  DURING  ITERATIVE  SOLUTION 

=  HOLDS  ORIGINAL  VALUE  OF  CBARK    DURING  ITERATIVE  SOLUTION 

=  HOLDS  ORIGINAL  VALUE  OF  CBARI  DURING  ITERATIVE  SOLUTION 

=  HOLDS  ORIGINAL  VALUE  OF  CBARJ  DURING  ITERATIVE  SOLUTION 

=  HOLDS  ORIGINAL  VALUE  OF  CBARK  DURING  ITERATIVE  SOLUTION 

=  HOLDS  ORIGINAL  VALUE  OF  CBARI  DURING  ITERATIVE  SOLUTION 

=  HOLDS  ORIGINAL  VALUE  OF  CBARJ  DURING  ITERATIVE  SOLUTION 

=  HOLDS  ORIGINAL  VALUE  OF  CBARK  DURING  ITERATIVE  SOLUTION 

=  VARIABLE  USED  IN  THE  ATHENB  GROUP  OF  SUBROUTINES  TO 
REPLACE  A  GROUP  OF  VARIABLES  FROM  BEING  RECALCULATED 
REPEATEDLY 

=  HOLDS  ORIGINAL  VALUE  OF  CI  DURING  ITERATIVE  SOLUTION 
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C  CJO        =  HOLDS  ORIGINAL  VALUE  OF  CJ  DURING  ITERATIVE  SOLUTION 

C 

C  CKO        =  HOLDS  ORIGINAL  VALUE  OF  CK  DURING  ITERATIVE  SOLUTION 

C 

C  OSCTST  =  VARIABLE  USED  TO  TEST  FOR  OSCILATORY  SOLUTION 

C 

C  RMSE      =  CALCULATED  ROOT  MEAN  SQUARE  ERROR 
C 

C  RMSEl     =  RMSE  VALUE  OF  PREVIOUS  ITERATION,   USED  TO  CHECK  AGAINST 

C  CONVERGENCE  FAILURE  DUE  TO  OSCILATORY  SOLUTIONS 

C 

C 

C 

C  SUBROUTINE  OUTPUT 

C 

C 

C  CT  =  TOTAL  NORMALITY,   RECALCULATED  AS  A  CHECK 

C 

C  CEC         =  CATION  EXCHANGE  CAPACITY,   RECALCULATED  AS  A  CHECK 

CALL  INPUT 
WRITE (4, 140) TITLE 
IDEBUG    =  0 
lUPPER    =  100 
ILOWER    =  -1 
JUPPER  =100 
JLOWER    =  -1 

CCTN  =  THETA*CTOTAL  +  BULKD*CECFIX 
EXPFRC  =  BULKD*CECFIX/CCTN 
SNPFRC  =  THETA*CTOTAL/CCTN 
N  =  MAXTAB 
RN  =  N 

DO  10  I  =  0,  N 
DO  9  J  =  0,N 
IF  ( I. LT. ILOWER)  GO  TO  9 
IFd.GT.IUPPER)  GO  TO  9 
IF (J. LT. JLOWER)   GO  TO  9 
IF (J. GT. JUPPER)  GO  TO  9 
NJ  =  N-J 
IF(I+J  -N)2,2,9 

2  CONTINUE 

TFRACJ  =  I/RN 
TFRACK  =  J/RN 

TFRACI  =  1.0  -  TFRACJ  -  TFRACK 
CALL  CCALC 
CALL  BLCKBX(I,J) 
CALL  FCALC 
IF(NJ-I)    3,    4,  3 

3  CONTINUE 

TABLE(1,NJ,I)   =  EPFRCK 
TABLE (1, I, N J)   =  EPFRCJ 
GO  TO  9 

4  CONTINUE 

TABLE (1, I, N J)   =  EPFRCJ 

9  CONTINUE 
10  CONTINUE 

DO  11  I  =  0,  N 

IF ( IDEBUG . LT . 1 ) WRITE (4,130)    (TABLE ( 1 , I , J) , J=0 , N) 
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11  CONTINUE 
100  FORMAT (7X, 12  (8X, 13, lOX)  ) 

110  FORMAT(2X,I2,2X,12(1X,F2.0,F5.2,F5.2,F5.2,3X) ) 

120  FORMAT (12X, •NA',3X, 'CA',3X, •MG',2(9X,  •NA',3X, •CA',3X, 'MG')) 

130  FORMAT (10G12. 5) 

140  FORMAT (A) 

STOP 

END 

C 
C 

C  ****************★******★*★***★**★*************^***^^^.**^^^^^^^^^^^^^^^ 

C 
C 
C 

SUBROUTINE  BLCKBX (III,  JJJ) 

C 
C 

C  **********************************************************^^^^^^^^^ 

C 
C 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI     IRF " 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKUO,  RKJK, 
&RKJKO,  RKKI,  RKKIO,  SNPFRC,  TABLE (4, 0 : 101, 0 : 101 ) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

1  CONTINUE 
ITERl  =  0 
ITER2  =  0 
CBARI 0  =  CBARI 
CBARJO  =  CBARJ 
CBARKO  =  CBARK 
5  CONTINUE 

CBARI  1  =  CBARI 
CBARJl  =  CBARJ 
CBARKl  =  CBARK 
IJ  =  0 
JK  =  0 
KI  =  0 
IJK  =  0 
10  CONTINUE 

IF(IJ.GE.l)   GO  TO  20 

IF (CCI.LE. 0.000001)   GO  TO  20 
IF(CCJ.LE. 0.000001)   GO  TO  20 
CALL  F 
CALL  GAMMA 
CALL  KCALC 
CALL  ITHENJ 
IF(IJ.GE.l)   GO  TO  20 

CALL  JTHENI 
IF(IJ.GE.l)   GO  TO  20 
GO  TO  30 
20  CONTINUE 

IF(IJK.GE.2)   GO  TO  90 
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30  CONTINUE 

IF(JK.GE.l)   GO  TO  40 

IF (CCJ.LE. 0.000001)  GO  TO  40 
IF (CCK.LE. 0.000001)   GO  TO  40 
CALL  F 
CALL  GAMMA 
CALL  KCALC 
CALL  JTHENK 
IF(JK.GE.l)   GO  TO  40 

CALL  KTHENJ 
IF(JK.GE.l)   GO  TO  40 
GO  TO  50 
40  CONTINUE 

IF(IJK.GE.2)   GO  TO  90 
50  CONTINUE 

IF(KI.EQ.l)  GO  TO  60 

IF(CCK.LE. 0.000001)   GO  TO  60 
IF (CCI.LE. 0.000001)   GO  TO  60 
CALL  F 
CALL  GAMMA 
CALL  KCALC 
CALL  KTHENI 
IF(KI.EQ.l)   GO  TO  60 
CALL  ITHENK 
60  CONTINUE 

IF(IJK-l)   100,   70,  90 
70  CONTINUE 

IF(CCI.LE. 0.000001)  GO  TO  130 
IF (CCJ.LE. 0.000001)   GO  TO  130 
IF (CCK.LE. 0.000001)   GO  TO  130 
IF(ITER2.GE.ITRLIM)   GO  TO  110 
ITER2  =  ITER2  +  1 
GO  TO  10 
90  CONTINUE 
ITER2  =  0 

RMSE  =  DSQRT ( ( (CBARI  -  CBARI1)**2  +   (CBARJ  -  CBARJ1)**2) 
&  +   (CBARK  -  CBARKl)**2)/3.0 
IF(RMSE.LE.TOL)   GO  TO  130 
C  IF(ITER1.GE.980) 

C        &WRITE(4,*)III,JJJ, •   CBARI        CBARI,   CBARJ,  CBARK 
IF(ITERl.GE.ITRLIM)   GO  TO  120 
ITERl  =  ITERl  +  1 
GO  TO  5 
100  CONTINUE 

IF(CCI.LE. 0.000001. AND. CCJ.LE. 0.000001)  GO  TO  130 
IF(CCI.LE. 0.000001. AND. CCK.LE. 0.000001)   GO  TO  130 
IF(CCJ.LE. 0.000001. AND. CCK.LE. 0.000001)   GO  TO  130 
IF(IDEBUG.GE.1)WRITE(4,*)   CCI,   CCJ,  CCK 
IF(IDEBUG.GE.1)WRITE(4,*)   CI,   CJ,  CK 
IF(IDEBUG.GE.1)WRITE(4,*)   CBARI, CBARJ,  CBARK 

WRITE (4,*) III,   JJJ, •     NONE  OF  THE  SIX  SUBROUTINES  CONVERGED' 

GO  TO  140 
110  CONTINUE 

WRITE (4,*) 'STABLE  FIRST  ESTIMATE  NOT  ACHIEVED' 

WRITE(4,*)III,  JJJ, 'AFTER',   ITER2,  '  ITERATIONS' 

GO  TO  140 
120  CONTINUE 

WRITE (4,*)    'STABLE  FINAL  ESTIMATE  NOT  ACHIEVED' 
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WRITE(4,*)  III,  JJJ, •  AFTER',  ITER2, '  ITERATIONS- 
GO  TO  140 

130  CONTINUE 

140  CONTINUE 

150  CONTINUE 
RETURN 
END 


C 
C 

C 
C 
C 

SUBROUTINE  CCALC 

C 
C 

c  *******************************************************^^*^^^^^^^^^ 

c 
c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI,  IRF 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKUO,  RKJK, 
&RKJKO,  RKKI,  RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 : 101) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

CCI  =  (CCTN*TFRACI)/ZI 

CCJ  =   (CCTN*TFRACJ) /Z J 

CCK  =   (CCTN*TFRACK) /ZK 

CI  =  CCI*SNPFRC/THETA 

CJ  =  CCJ*SNPFRC/THETA 

CK  =  CCK*SNPFRC/THETA 

CBARI  =  CCI*EXPFRC/BULKD 

CBARJ  =  CCJ*EXPFRC/BULKD 

CBARK  =  CCK*EXPFRC/BULKD 

CALL  DNMNTR 

RETURN 
END 


C 

C 

c 
c 
c 

SUBROUTINE  DNMNTR 

C 

c 

c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ  CCK 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
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&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI,  IRF, 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,  RKKI,  RKKIO,   SNPFRC,   TABLE (4 , 0 : 101, 0 : 101 ) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

IF(IRF-2)  1,3,2 

1  CONTINUE 

DENOM  =  CBARI  +  CBARJ  +  CBARK 
GO  TO  7 

2  CONTINUE 

DENOM  =   ((ZI  +  1.0)*CBARI  +    (ZJ  +  1.0)*CBARJ  + 
&  (ZK  +  1.0)*CBAARK)/2.0 

GO  TO  7 

3  CONTINUE 

CEC  =  ZI*CBARI  +  ZJ*CBARJ  +  ZK*CBARK 
IF(ISBMKR.EQ.l)  GO  TO  4 
IF(ISBMKR.EQ.2)  GO  TO  6 
IF(ISBMKR.EQ.3)  GO  TO  4 
IFdSBMKR.EQ.  4)  GO  TO  5 
IF(ISBMKR.EQ.5)  GO  TO  6 
IFdSBMKR.EQ.  6)   GO  TO  5 

4  CONTINUE 

DENOM  =  CEC*((ZJ**ZI)*(ZI**(-ZJ)))**(1.0/(ZJ-ZI)) 
GO  TO  7 

5  CONTINUE 

DENOM  =  CEC* ( (ZK**Z J) * (Z J** (-ZK) ) ) ** (1 . 0/ (ZK-ZJ) ) 
GO  TO  7 

6  CONTINUE 

DENOM  =  CEC* ( (ZI**ZK) * (ZK** (-ZI) ) ) ** (1 . 0/ (ZI-ZK) ) 
GO  TO  7 
7  CONTINUE 
RETURN 
END 


C 

C 

c  *****************************************************^^^^^^^^^^^^^^ 
C 

c 
c 

SUBROUTINE  F 

C 

C 

c  **************************************************^*^^^^^^^^^^^^^^^ 

C 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

CQMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ  CCK 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI  IRF 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK 
&RKJKO,   RKKI,   RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 • 101 ) 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ  ZK 


FI  =  1.0 
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FJ  =  1.0 
FK  =  1.0 
RETURN 
END 

C 
C 

c  ******************************************************************* 

c 
c 
c 

SUBROUTINE  FCALC 

C 

c 

c   ************ ********************************************.*********** 

c 
c 

c 

IMPLICIT  DOUBLE  PRECISION  {A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK,  IDAVIS,   IDEBUG,    IJ,    IJK,   KI,  IRF, 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,   RKKI,   RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 : 101) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,    ZI,   ZJ,  ZK 

TEQUIV  =  ZI*CBARI  +  ZJ*CBARJ  +  ZK*CBARK 

EPFRCI  =   (ZI*CBARI) /TEQUIV 

EPFRCJ  =   (ZJ*CBARJ) /TEQUIV 

EPFRCK  =   (ZK*CBARK) /TEQUIV 

RETURN 

END 

C 

c 

c  *********************************************************^^^^^^^^^^ 

c 
c 
c 

SUBROUTINE  GAMMA 

C 

c 

c  *********************************************^,*****^^^^^^^^^^^^^^^^ 

c 
c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK. 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,  EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,    GAMMAJ,    GAMMAK,  IDAVIS,    IDEBUG,    I  J,    IJK,   KI  IRF 
&ISBMKR,    ITRLIM,    JK,     MAXTAB,   RKIJ,    RKIJO,  RKJK, 
&RKJKO,   RKKI,   RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0  : 101 ) 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

IF (IDAVIS)  10,10,5 
5  CONTINUE 

ZISQR  =  ZI**2.0 
ZJSQR  =  ZJ**2.0 
ZKSQR  =  ZK**2,0 
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SI  =   (CI*ZISQR  +  CJ*ZJSQR  +  CK*ZKSQR) /2 . 0 
RTSI  =  DSQRT(SI) 

DVSBSS  =  -   ((RTSI/(1.0  +  RTSI) ) -0 . 2*SI) /2 .  0 
GAMMAI  =  DEXP(ZISQR*DVSBSS) 
GAMMAJ  =  DEXP (ZJSQR* DVSBSS) 
GAMMAK  =  DEXP (ZKSQR* DVSBSS) 
10  CONTINUE 

GAMMAI  =1,0 
GAMMAJ  =1.0 
GAMMAK  =  1.0 
RETURN 
END 


C     **************  ************* 


*************************** 


SUBROUTINE  INPUT 


*****************************^^^ 


C 
C 

c 
c 

c  *********************************** 

c 
c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   I J,  IJK, 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,   RKKI,   RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 : 101 )  , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,  ZJ, 

READ (10, 1000)  TITLE 


KI,  IRF, 


ZK 


READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1010) 
READ (10, 1020) 
READ (10, 1020) 
READ (10, 1020) 
READ (10, 1020) 
READ (10, 1020) 
WRITE (*, 1000) 
WRITE (*, 1010) 
WRITE (*,  1010) 
WRITE (*,  1010) 
WRITE (*,  1010) 
WRITE (*,  1010) 


BULKD 

CECFIX 

CTOTAL 

THETA 

ZI 

ZJ 

ZK 

RKIJO 
RKJKO 
RKKIO 
TOL 

ITRLIM 
MAXTAB 
IRF 

IDAVIS 

IDEBUG 

TITLE 

BULKD 

CECFIX 

CTOTAL 

THETA 

ZI 


162 


WRITE (*, 1010)  ZJ 

WRITE (*, 1010)  ZK 

WRITE (*, 1010)  RKIJO 

WRITE (*, 1010)  RKJKO 

WRITE (*, 1010)  RKKIO 

WRITE (*, 1010)  TOL 

WRITE (*, 1020)  ITRLIM 

WRITE (*, 1020)  MAXTAB 

WRITE (*, 1020)  IRF 

WRITE (*, 1020)  IDAVIS 

WRITE (*, 1020)  IDEBUG 
1000  FORMAT (A) 
1010  FORMAT (G15,  7) 
1020  FORMAT (110) 
1030  FORMAT (7X, A) 

RETURN 

END 

C 

c 

c    *************** ****************************** *********************** 

c 

c 

c 

SUBROUTINE  ITHENJ 

C 

C 

C  *********************************************************v^*^k.*j^^*^^^^ 

c 
c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,    IJ,   IJK,   KI,  IRF, 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,  RKKI,  RKKIO,  SNPFRC,  TABLE (4, 0 : 101, 0 : 101) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

CBARI  0  =  CBARI 

CBARJO  =  CBARJ 

CBARI 1  =  CBARI 

CBARJl  =  CBARJ 

CIO  =  CI 

CJO  =  CJ 

ITER  =  0 

RMSE  =0.01 
10  CONTINUE 

ITER  =  ITER  +  1 

RMSEl  =  RMSE 

CBARI2  =  CBARIl 

CBARJ2  =  CBARJl 

CBARIl  =  CBARI 
CBARJl  =  CBARJ 

CBICI  =   ((((CBARJ/CJ)**ZI)*(DEN0M**(ZJ-ZI)))/RKIJ)**(1  0/ZJ) 

CI  =  CCI/ (BULKD*CBICI  +  THETA) 

CBARI  =  CCI/((THETA/CBICI)   +  BULKD) 

CJ  =   (CTOTAL  -  ZI*CI  -  ZK*CK) /ZJ 

CBARJ  =   (CECFIX  -  ZI*CBARI  -  ZK*CBARK) /Z J 
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IF(IDEBUG.GE.l)   WRITE (4, *) ' ITHENJ' 
IF(IDEBUG.GE.l)  WRITE (4,*)  ITER 

IF(IDEBUG.GE.l)  WRITE (4, *) ' CI,   CJ,   CK     ',CI,   CJ,  CK 
IF(IDEBUG.GE.l)  WRITE  (4,  *)  'CBARI,   CBARJ,   CBARK  SCBARI,  CBARJ, 

CBARK 

CALL  DNMNTR 

IF{CJ)   60,   60,  20 
20  CONTINUE 

IF (CBARJ)   60,   60,  30 
30  CONTINUE 

RMSE  =  DSQRT (( (CBARI  -  CBARI1)**2  +    (CBARJ  -  CBARJl) **2) /2) 

IF(RMSE  -  TOL)   50,   50,  40 

40  CONTINUE 
IF(ITER-2)    49,   41,  41 

41  CONTINUE 

OSCTST  =  DSQRT (( (CBARI -CBARI2)**2  +   (CBARJ-CBARJ2) **2) /2) 
IF(OSCTST  -  TOL)  43,43,42 

42  CONTINUE 
IF(RMSEl-RMSE)  43,49,49 

43  CONTINUE 

CBARJ  =   (CBARJ  +  CBARJl )/2 

CJ  =   (CCJ-(CBARJ*BULKD) )/THETA 

CBARI  =   (CECFIX  -  ZJ*CBARJ  -  ZK*CBARK) /ZI 

CI  =(CTOTAL  -  ZJ*CJ  -  ZK*CK)/ZI 

RMSE  =0.01 

49  CONTINUE 

IFdTRLIM  -  ITER)    60,   10,  10 

50  CONTINUE 
IJ  =  1 

GO  TO  70 
60  CONTINUE 

CI  =  CIO 
CJ  =  CJO 
CBARI  =  CBARI  0 
CBARJ  =  CBARJO 
CALL  DNMNTR 
70  CONTINUE 

IJK  =  IJ  +  JK  +  KI 
RETURN 
END 

C 

c 

c 
c 
c 

SUBROUTINE  ITHENK 

C 
C 

C 

c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI  IRF 
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&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,  RKKI,  RKKIO,   SNPFRC,  TABLE (4, 0 : 101, 0 : 101) , 
STFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,    ZJ,  ZK 

CBARIO  =  CBARI 

CBARKO  =  CBARK 

CBARI 1  =  CBARI 

CBARKl  =  CBARK 

CIO  =  CI 

CKO  =  CK 

ITER  =  0 

RMSE  =0.01 
10  CONTINUE 

ITER  =  ITER  +  1 

RMSEl  =RMSE 

CBARI 2  =  CBARI 1 

CBARK2  =  CBARKl 

CBARI  1  =  CBARI 

CBARKl  =  CBARK 

CBICI  =   (RKKI*( (CBARK/CK)**ZI)*(DENa^**(ZK-ZI) ) )**(1  0/ZK) 

CI  =  CCI/ (BULKD*CBICI  +  THETA) 

CBARI  =  CCI/ ( (THETA/CBICI)   +  BULKD) 

CK  =   (CTOTAL  -  ZI*CI  -  Z J*CJ) /ZK 

CBARK  =   (CECFIX  -  ZI*CBARI  -  ZJ*CBARJ) /ZK 

CALL  DNMNTR 

IF(IDEBUG.GE.1)WRITE(4,*)  'ITHENK' 
IF(IDEBUG.GE.1)WRITE(4,*)  ITER 

IF(IDEBUG.GE.1)WRITE(4,*) -CI,   CJ,   CK     •,CI,   CJ,  CK 
IF(IDEBUG.GE.1)WRITE(4,*) -CBARI,   CBARJ,   CBARK  CBARI, 
&CBARJ,  CBARK 
IF(CK)   60,   60,  20 
20  CONTINUE 

IF (CBARK)   60,   60,  30 
30  CONTINUE 

RMSE  =  DSQRT (( (CBARI  -  CBARI1)**2  +  (CBARK  -  CBARKl) **2) /2 
IF  (RMSE  -  TOL)    50,   50,  40 

40  CONTINUE 
IF(ITER-2)   49,   41,  41 

41  CONTINUE 

OSCTST  =  DSQRT (( (CBARI-CBARI2)**2  +   (CBARK-CBARK2) **2) 
IF(OSCTST  -  TOL)  43,43,42 

42  CONTINUE 
IF(RMSEl-RMSE)  43,49,49 

43  CONTINUE 

CBARK  =   (CBARK  +  CBARKl) /2 

CK  =   (CCK-(CBARK*BULKD) ) /THETA 

CBARI  =   (CECFIX  -  ZJ*CBARJ  -  ZK*CBARK) /ZI 

CI  = (CTOTAL  -  ZJ*CJ  -  ZK*CK)/ZI 

RMSE  =  0.01 

49  CONTINUE 

IF (ITRLIM  -  ITER)    60,   10,  10 

50  CONTINUE 
KI  =  1 

GO  TO  70 
60  CONTINUE 
CI  =  CIO 
CK  =  CKO 
CBARI  =  CBARIO 
CBARK  =  CBARKO 
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CALL  DNMNTR 
70  CONTINUE 

IJK  =  IJ  +  JK  +  KI 
RETURN 
END 

C 
C 

c 
c 
c 

SUBROUTINE  JTHENI 

C 
C 

C  ******************************************** 

c 
c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK,  IDAVIS,   IDEBUG,    IJ,    IJK,   KI,  IRF 
&ISBMKR,    ITRLIM,   JK,     MAXTAB,   RKIJ,   RKUO,  RKJK, 
&RKJKO,   RKKI,  RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 : 101) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

CBARI 0  =  CBARI 

CBARJO  =  CBARJ 

CBARI 1  =  CBARI 

CBARJl  =  CBARJ 

CIO  =  CI 

CJO  =  CJ 

ITER  =  0 

RMSE  =0.01 
10  CONTINUE 

ITER  =  ITER  +  1 

RMSEl  =  RMSE 

CBARI2  =  CBARIl 

CBARJ2  =  CBARJl 

CBARIl  =  CBARI 

CBARJl  =  CBARJ 

CBJCJ  =  (RKIJ* ((CBARI/CI)**ZJ)*( (DENOM) **(ZI-ZJ)))**(1  0/ZI) 
CJ  =  CCJ/ (BULKD*CBJCJ  +  THETA)  ^  -u/^x; 

CBARJ  =  CCJ/ ( (THETA/CBJCJ)   +  BULKD) 

CI  =   (CTOTAL  -  ZJ*CJ  -  ZK*CK)/ZI 

CBARI  =   (CECFIX  -  ZJ*CBARJ  -  ZK*CBARK) /ZI 

CALL  DNMNTR 

IF(IDEBUG.GE.1)WRITE{4,*)  'JTHENI' 
IF(IDEBUG.GE.1)WRITE(4,*)  ITER 

IF(IDEBUG.GE.1)WRITE(4,*) 'CI,   CJ,   CK     ',CI,   CJ,  CK 
IF(IDEBUG.GE.1)WRITE(4,*) 'CBARI,   CBARJ,   CBARK  ' 
& CBARI,   CBARJ,  CBARK 
IF (CI)    60,   60,  20 
20  CONTINUE 

IF (CBARI)   60,   60,  30 
30  CONTINUE 

RMSE  =  DSQRT  ((  (CBARI  -  CBARIl)  **2  +  (CBARJ  -  CBARJl)  **2) /2) 
IF (RMSE  -  TOL)   50,   50,   40  '  i/^> 
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40  CONTINUE 
IF(ITER-2)   49,   41,  41 

41  CONTINUE 

OSCTST  =  DSQRT( ( (CBARI-CBARI2)**2  +   (CBARJ-CBARJ2) **2) /2) 
IF(OSCTST  -  TOL)  43,43,42 

42  CONTINUE 
IF(PMSEl-RMSE)  43,49,49 

43  CONTINUE 

CBARI  =   (CBARI  +  CBARIl)/2 

CI  =   (CCI-(CBARI*BULKD) )/THETA 

CBARJ  =   (CECFIX  -  ZI*CBARI  -  ZK*CBARK) /ZJ 

CJ  ={CTOTAL  -  ZI*CI  -  ZK*CK)/ZJ 

RMSE  =0.01 

49  CONTINUE 

IFdTKLIM  -  ITER)    60,   10,  10 

50  CONTINUE 
IJ  =  1 

GO  TO  70 
60  CONTINUE 
CI  =  CIO 

CJ  =  CJO 

CBARI  =  CBARI 0 

CBARJ  =  CBARJO 

CALL  DNMNTR 
70  CONTINUE 

IJK  =  IJ  +  JK  +  KI 

RETURN 

END 

C 
C 

C 
C 
C 

SUBROUTINE  JTHENK 

C 
C 

C 

c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI,  IRF 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,   RKKI,   RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 : 101 ) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

CBARJO  =  CBARJ 

CBARKO  =  CBARK 

CBARJl  =  CBARJ 

CBARKl  =  CBARK 

CJO  =  CJ 

CKO  =  CK 

ITER  =  0 

RMSE  =0.01 
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10  CONTINUE 

ITER  =  ITER  +  1 
RMSEl  =  RMSE 
CBARJ2  =  CBARJl 
CBARK2  =  CBARKl 
CBARJl  =  CBARJ 
CBARKl  =  CBARK 

CBJCJ  =   ( ( ( (CBARK/CK) **ZJ) * ( (DENC^) 
&** (ZK-Z J) ) ) /RKJK) ** (1 . 0/ZK) 
CJ  =  CCJ/ (BULKD*CBJCJ  +  THETA) 
CBARJ  =  CCJ/ ( (THETA/CBJCJ)   +  BULKD) 
CK  =   (CTOTAL  -  ZI*CI  -  ZJ*CJ) /ZK 
CBARK  =   (CECFIX  -  ZI*CBARI  -  ZJ*CBARJ) /ZK 
CALL  DNMNTR 

RMSE  =  DSQRT {( (CBARJ  -  CBARJl) **2  +    (CBARK  -  CBARKl) **2) /2) 
IF(IDEBUG.GE.l)  WRITE (4,*)  'JTHENK' 
IF(IDEBUG.GE,1)WRITE(4,*)  ITER 

IF(IDEBUG.GE.1)WRITE(4,*) -CI,   CJ,   CK     ',CI,   CJ,  CK 
IF(IDEBUG.GE.1)WRITE(4,*) 'CBARI,   CBARJ,   CBARK  ', 
&CBARI,   CBARJ,  CBARK 
IF(CK)   60,   60,  20 
20  CONTINUE 

IF  (CBARK)   60,   60,  30 
30  CONTINUE 

RMSE  =  DSQRT {( (CBARJ  -  CBARJl) **2  +    (CBARK  -  CBARKl) **2) /2) 
IF (RMSE  -  TOL)   50,   50,  40 

40  CONTINUE 
IF(ITER-2)   49,   41,  41 

41  CONTINUE 

OSCTST  =  DSQRT (((CBARJ-CBARJ2)**2  +    (CBARK-CBARK2) **2) /2) 
IF(OSCTST  -  TOL)  43,43,42 

42  CONTINUE 
IF(RMSEl-RMSE)  43,49,49 

43  CONTINUE 

CBARK  =   (CBARK  +  CBARKl) /2 

CK  =   (CCK-(CBARK*BULKD) ) /THETA 

CBARJ  =   (CECFIX  -  ZI*CBARI  -  ZK*CBARK) /ZJ 

CJ  = (CTOTAL  -  ZI*CI  -  ZK*CK)/ZJ 

RMSE  =0,01 

49  CONTINUE 

IF(ITRLIM  -  ITER)    60,   10,  10 

50  CONTINUE 
JK  =  1 

GO  TO  70 
60  CONTINUE 
CJ  =  CJO 

CK  =  CKO 

CBARJ  =  CBARJO 

CBARK  =  CBARKO 

CALL  DNMNTR 
70  CONTINUE 

IJK  =  IJ  +  JK  +  KI 

RETURN 

END 

C 
C 
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C  ******************************************************************* 

c 
c 
c 

SUBROUTINE  KCALC 

C 

c 

C 

c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI,  IRF, 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,   RKKI,  RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 : 101) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

RKIJ  =  RKIJO  *  ((  FI/GAMMAI  ) **Z J  )  * 
&((  GAMMAJ/FJ  )**ZI  ) 

RKJK  =  RKJKO  *  ( (  FJ/GAMMAJ  ) **ZK  )  * 
& ( (  GAMMAK/FK  ) **Z J  ) 

RKKI  =  RKKIO  *  ((  FK/GAMMAK  ) **ZI  )  * 
&  ( (  GAt^I/FI  )  **ZK  ) 

RETURN 

END 

C 

C 

c  ******************************************************************* 

C 

c 

C 

SUBROUTINE  KTHENI 

C 
C 

C  ***********************************************************^^**^^^^^^ 
C 

c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI,  IRF 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,   RKKI,   RKKIO,   SNPFRC,   TABLE (4, 0 : 101, 0 : 101) 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

CBARIO  =  CBARI 

CBARKO  =  CBARK 

CBARI  1  =  CBARI 

CBARKl  =  CBARK 

CIO  =  CI 

CKO  =  CK 

ITER  =  0 

RMSE  =0.01 
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10  CONTINUE 

RMSEl  =  RMSE 
ITER  =  ITER  +  1 
CBARI2  =  CBARIl 
CBARK2  =  CBARKl 
CBARIl  =  CBARI 
CBARKl  =  CBARK 
CALL  DNMNTR 

CBKCK  =  ((((CBARI/CI)**ZK)*(DENOM**(ZI-ZK)))/RKKI)**(1.0/ZI) 

CK  =  CCK/ (BULKD*CBKCK  +  THETA) 

CBARK  =  CCK/ ( (THETA/CBKCK)   +  BULKD) 

CI  =   (CTOTAL  -  ZJ*CJ  -  ZK*CK)/ZI 

CBARI  =   (CECFIX  -  ZJ*CBARJ  -  ZK*CBARK) /ZI 

CALL  DNMNTR 

IF(IDEBUG.GE.1)WRITE(4,*)  'KTHENI' 

IF(IDEBUG.GE.1)WRITE(4,*) 'CI,   CJ,   CK     ' ,CI,   CJ,  CK 
IF(IDEBUG.GE.1)WRITE(4,*) 'CBARI,   CBARJ,   CBARK  ', 
& CBARI,   CBARJ,  CBARK 
IF (CI)   60,   60,  20 

20  CONTINUE 

IF (CBARI)   60,   60,  30 

30  CONTINUE 

RMSE  =  DSQRT (( (CBARI  -  CBARIl) **2  +   (CBARK  -  CBARKl) **2) /2) 
IF(IDEBUG.GE.1)WRITE(4,*)    ' KTHENI ', ITER, RMSE,   RMSEl, CBARI 
IF  (IDEBUG.GE.l)WRITE (4, *)  ITER 

IF(IDEBUG.GE.1)WRITE(4,*) 'CI,   CJ,   CK     ',CI,   CJ,  CK 
IF (IDEBUG.GE.l) WRITE (4,*) 'CBARI,   CBARJ,   CBARK  ', 
& CBARI,   CBARJ,  CBARK 
IF  (RMSE  -  TOL)   50,   50,  40 

40  CONTINUE 
IF(ITER-2)   49,   41,  41 

41  CONTINUE 

OSCTST  =  DSQRT (( (CBARI-CBARI2)**2  +   (CBARK-CBARK2) **2) /2) 
IF(OSCTST  -  TOL)  43,43,42 

42  CONTINUE 
IF(RMSEl-RMSE)  43,49,49 

43  CONTINUE 

CBARI  =   (CBARI  +  CBARIl) /2 

CI  =   (CCI-(CBARI*BULKD) ) /THETA 

CBARK  =   (CECFIX  -  ZI*CBARI  -  Z J*CBARJ) /ZK 

CK  = (CTOTAL  -  ZI*CI  -  Z J*CJ) /ZK 

RMSE  =0.01 

49  CONTINUE 

IF(ITRLIM  -  ITER)    60,   10,  10 

50  CONTINUE 
KI  =  1 

GO  TO  70 
60  CONTINUE 
CI  =  CIO 

CK  =  CKO 

CBARI  =  CBARI 0 

CBARK  =  CBARKO 

CALL  DNMNTR 
70  CONTINUE 

IJK  =  IJ  +  JK  +  KI 

RETURN 

END 
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C 

c 

c 
c 
c 

SUBROUTINE  KTHENJ 

C 
C 

Q  ***************************************** 

c 
c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COIMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI,  IRF, 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,  RKKI,  RKKIO,   SNPFRC,  TABLE (4, 0 : 101, 0 : 101) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 
CBARI 0  =  CBARI 

CBARJO  =  CBARJ 

CBARKO  =  CBARK 

CBARJl  =  CBARJ 

CBARKl  =  CBARK 

CJO  =  CJ 

CKO  =  CK 

ITER  =  0 

RMSE  =0.01 
10  CONTINUE 

ITER  =  ITER  +  1 

RMSEl  =  RMSE 

CBARJ2  =  CBARJl 

CBARK2  =  CBARKl 

CBARJl  =  CBARJ 

CBARKl  =  CBARK 

CBKCK  =   (RKJK* ( (CBARJ/CJ) **ZK) * (DENOM** (ZJ-ZK) ) ) ** (1 . 0/ZJ) 

CK  =  CCK/ (BULKD*CBKCK  +  THETA) 

CBARK  =  CCK/ ( (THETA/CBKCK)   +  BULKD) 

CJ  =   (CTOTAL  -  ZI*CI  -  ZK*CK)/ZJ 

CBARJ  =   (CECFIX  -  ZI*CBARI  -  ZK*CBARK) /Z J 

CALL  DNMNTR 

RMSE  =  DSQRT (( (CBARJ  -  CBARJl) **2  +   (CBARK  -  CBARKl) **2) /2) 
IF(IDEBUG,GE.1)WRITE{4,*)  'KTHENJ' 
IF(IDEBUG.GE.1)WRITE(4,*)  ITER 

IF(IDEBUG.GE.1)WRITE(4,*) 'CI,   CJ,   CK     ',01,   CJ,  CK 
IF{IDEBUG.GE.1)WRITE(4,*) 'CBARI,   CBARJ,   CBARK  ', 
&CBARI,    CBARJ,  CBARK 
IF(CJ)    60,   60,  20 
20  CONTINUE 

IF (CBARJ)   60,   60,  30 
30  CONTINUE 

RMSE  =  DSQRT (( (CBARJ  -  CBARJl) **2  +  (CBARK  -  CBARKl ) **2) /2) 
IF (RMSE  -  TOL)   50,   50,  40 

40  CONTINUE 
IF(ITER-2)   49,   41,  41 

41  CONTINUE 
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OSCTST  =  DSQRT( ( (CBARJ-CBARJ2)**2  +   (CBARK-CBARK2) **2) /2) 
IF(OSCTST  -  TOL)  43,43,42 

42  CONTINUE 
IF(RMSEl-RMSE)  43,49,49 

43  CONTINUE 

CBARJ  =   (CBARJ  +  CBARJl)/2 

CJ  =   (CCJ-(CBARJ*BULKD) )/THETA 

CBARK  =   (CECFIX  -  ZI*CBARI  -  Z J*CBARJ) /ZK 

CK  =(CTOTAL  -  ZI*CI  -  ZJ*CJ) /ZK 

RMSE  =0.01 

49  CONTINUE 

IFdTRLIM  -  ITER)   60,  10,  10 

50  CONTINUE 
JK  =  1 

GO  TO  70 
60  CONTINUE 
CJ  =  CJO 

CK  =  CKO 

CBARJ  =  CBARJO 

CBARK  =  CBARKO 

CALL  DNMNTR 
70  CONTINUE 

IJK  =  IJ  +  JK  +  KI 

RETURN 

END 

C 
C 

C  *******************************************************^*^^^^^^^^^^ 

C 

C 

c 

SUBROUTINE  KWRITE 

C 

c 

c 
c 
c 
c 
c 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENQM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI  IRF 
&ISBMKR,    ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJK0,   RKKI,   RKKIO,   SNPFRC,   TABLE (4, 0 : 101 , 0 • 101) 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 

TOT  =  CBARI  +  CBARJ  +  CBARK 

TN  =  ZI*CI  +  ZJ*CJ  +  ZK*CK 

CEC=  ZI*CBARI  +  ZJ*CBARJ  +  ZK*CBARK 

CBARI  =  CBARI  +  O.ld-14 

CBARJ  =  CBARJ  +  O.ld-14 

CBARK  =  CBARK  +  0,ld-14 

CI  =  CI  +  O.ld-14 

CJ  =  CJ  +  O.ld-14 

CK  =  CK  +  O.ld-14 

RIJ  =   ( (CI/CBARI) **ZJ) * ( (CBARJ/CJ) **ZI) * (TOT** (ZJ-ZI) ) 
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RJK  =   ( (CJ/CBARJ) **ZK) * ( (CBARK/CK) **ZJ)  * (TOT**  (ZK-ZJ)  ) 
RKI  =  ((CK/CBARK)**ZI)*((CBARI/CI)**ZK)*(TOT**(ZI-ZK)) 
WRITE{4,100)TFRACI,   TFRACJ,   TFRACK,   RIJ,   RJK,  RKI 
100    FORMAT (IX,  6F1 2, 5) 
RETURN 
END 


C 
C 

C  ******************************************************************** 

c 
c 
c 

SUBROUTINE  OUTPUT 

C 
C 

Q  *************************************************.****:**********^^t*^ 

C 
C 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STEVES/BULKD,   CBARI,   CBARJ,   CBARK,   CCI,   CCJ,  CCK, 
&CCTN,   CECFIX,   CI,   CJ,   CK,   CTOTAL,  DENOM, 
&FI,   FJ,   FK,   EPFRCI,   EPFRCJ,   EPFRCK,  EXPFRC, 
&GAMMAI,   GAMMAJ,   GAMMAK, IDAVIS,   IDEBUG,   IJ,   IJK,   KI,  IRF, 
&ISBMKR,   ITRLIM,   JK,     MAXTAB,   RKIJ,   RKIJO,  RKJK, 
&RKJKO,  RKKI,  RKKIO,   SNPFRC,  TABLE (4, 0 : 101, 0 : IQl) , 
&TFRACI,   TFRACJ,   TFRACK,   THETA,   TITLE,   TOL,   ZI,   ZJ,  ZK 


CT  =  ZI*CI  +  ZJ*CJ  +ZK*CK 

CEC  =  ZI*CBARI  +  ZJ*CBARJ  +  ZK*CBARK 

CCI  =  THETA*CI  +  BULKD*CBARI 

CCJ  =  THETA*CJ  +  BULKD*CBARJ 

CCK  =  THETA*CK  +  BULKD*CBARK 

RETURN 

END 


APPENDIX  C 

OPERATING  PROGRAM  BLCKBX 
Program  BLCKBX  was  written  in  FORTRAN  77  and  implemented 
on  a  Digital  Equipment  Corporation  VAX  11/750  minicomputer. 
To  run  the  program,   an  input  file  with  physical,  chemical, 
and  computational  parameters  is  written.  Below  is  a  listing 
of  an  input  file  which  would  instruct  Program  BLCKBX  to 
created  a  50x50x50  ternary  exchange  isotherm  using  the  mole 
fraction  exchanger  phase  reference  function. 


MOLE  FRACTION 
1.161239 
0. 4494000E-02 
0.1055000E-01 
0.2510000 
1.000000 
2.000000 
2.000000 
7.215056 
0.4120000E-02 

2.159399 
O.lOOOOOOE-08 
100000 
50 
1 
0 
0 


REFERENCE  FUNCTION 

Bulk  Density  (kg/dec3) 
Cation  Exchange  Capacity  (mol(+)/kg) 
Total  Normality  (mol(+)/dec3) 
Moisture  content  (dec3/dec3) 
Valence  for  the  1st  ion 
Valence  for  the  2nd  ion 
Valence  for  the  3rd  ion 

Thermodynamic  exchange  constant  for  KEX12 
Thermodynamic  exchange  constant  for  KEX23 
Thermodynamic  exchange  constant  for  KEX31 
Tolerance  for  root  mean  square  error 
Iterative  limit 
Maximvun  index  in  table  {0..x) 

Exchange  Phase  Reference  Function=Mole  Fraction 
Implement  Davies  Equation 
Print  Debugging  Write  Statements 


When  running  the  program  interactively,   the  program  must 
first  be  compiled  using  the  commands 

$  FORTRAN  BLPCBOX  <return> 
$  LINK  BLKBOX  <return> 
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If  the  file  BLKBX.exe  is  maintained  in  the  user's 
directory,   the  two  lines  above  do  not  need  to  be  repeated. 
Then  the  program  is  the  begun  by  typing 
$  RUN  BLKBOX  <return> 
The  program  then  requests 

INPUT  FILE  NAME 

The  user  then  types 
input  file  name  <return> 
The  computer  then  inquires 
OUTPUT  FILE  NAME 

and  the  user  designates  his  or  her  choice 
output  file  name  <return> 
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